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ABSTRACT 

Under the NASA Contract NAS8-38609.D.O.6, which is the first year effort of the 
current model improvement efforts for predicting spary combustion processes in 
liquid-fueled rockets, four major tasks including atomization models, PDF droplet 
dispersion models, coalescence and breakup models, as well as dense-spray turbu- 
lence modulation effects have been incorporated into the CFD code MAST ( Multi- 
phase All-Speed Transient program ). 

The atomization model implemented is the “blob injection” model of Reitz in- 
volving secondary breakup mechanism. Two breakup models using the Taylor Anal- 
ogy Breakup ( TAB ) concept and a wave instability concept were compared with si- 
multaneous incorporation of an existing stochastic collision-coalescence model. Two 
dense modulation models based on the continuum approach of Chen and Fashola, 
and the stochastic approach of the Los Alamos group were also implemented and 
compared. To improve the computational efficiency, a parcel probability dense func- 
tion PDF) tracking method accounting for the dispersion within the numerical par- 
ticles was improved and implemented. 

Detailed formulations as well as validation studies are described in this report. 
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II 


Governing Equations and 
Physical Models 

II. 1 Introduction 

In this Chapter, all governing equations used to describe the spray combustions are 
presented. The gas-phase and liquid-phase processes are modeled by a system of 
unsteady, two-dimensional (axisymmetric) equations. The gas-phase equations are 
written in an Eulerian coordinate whereas the liquid-phase is presented in Lagrangian 
coordinates. The two-way coupling between the two phases is described by the 
interaction source terms which represent the rates of momentum, mass and heat 
exchange. A two-equation k - t turbulence model with particle modulation effect 
is used to characterize the gas phase flow properties. Particle motion equation 
is the simplified B-B-0 equation. Poly dispersed particle dynamics and particle 
turbulent dispersion are modelled using a Mote Carlo method. Parcel PDF model is 
used to improve computational efficiency. Droplet evaporation and heat transfer are 
calculated using Frossling correlation and Ranz-Marshall correlation respectively. A 
turbulent diffusion flame and single step chemical reaction model is used for spray 
combustion. Dense spray effects are accounted for by droplet breakup and collision 
models. 
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II. 2 Gas Phase Equations 

11*2.1 Mean Flow Equation 

The density- weighted conservation equations of mass, momentum, energy and mass 
fraction variables in an Eulerian coordinate can be written as followings: 


dp d 

dt + d7^ pUt ^ = Sm ’ 1 


(III) 


dpU, 9 dP d Tii n 

— 0 - I- S Ui + 5 Uli / 


dt dxj 


dx{ dxj 


( 112 ) 


dpHt ,d dP dQ j „ 

dt + dxj ( pUjHt ^ ~ dt ~7h~ + Sk + ShJ + RfuHc 


(11.3) 


dpY f d , dN f 


(11.4) 


dpK 

dt 

In the above equations, 


+ ^p v ’ Y * ) = ~ihf ~ s (H.5) 

p is the ensemble averaged density of the mixture, U t is 


the i component of the density- weighted (Favre) mean velocity, H t = H + \UiU t is 
the density-weighted mean gas total enthalpy and H is the gas static enthalpy, Y } 
and Y 02 are the density-weighted mean mass fractions of fuel and oxidizer, and P is 
the mean pressure. S m ,i, S u>: i and Sh,i represent interaction terms (or source terms) 
of spray to be defined later. R fu is the combustion rate in eddy breakup model, H com 
is the combustion heat of fuel vapor and s is stoichiometric coefficient for oxidizer 
and fuel. S h is the energy dissipation term and diffusion enthalpy is neglected in 


energy equation for multi-species flow, corresponding to unit Lewis number. 


4 


To account for turbulent flows, an eddy viscosity type model was employed. Thus 


, aUi du j. 

(II.6) 

^ Kp r t Prjdxj 

(11-1) 

a JM 8Yj 
NlJ - “ ( Sc + Sc>dx, 

(II.8) 

JV . = _(■£.+ '" ) dY » 
Sc Sc, 1 dx, 

(II.9) 

- 2. . . a ,3U, 

Sul ~ ^ dx^ dxj 

(IUO) 


S h = 


+ 


_d_ 

dxj 

_d_ 

dxj 




i 

Pr 

dUj 


dxj 


( n ' JWj ' i ^ ® t tj ®H±\ 
ip,U ' dx, 1 dxj 


( 11 . 11 ) 


in which fi t is described by the k — e two-equation model. The last two terms in Sh 
are neglected conveniently with little error. 


II. 2. 2 Equation of State 


In this study, the gas phase is assumed to be the mixture of ideal gases. The equation 
of state for the mixture is 


P = pRT 


( 11 . 12 ) 
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The gas constant R is related to the universal gas constant R con (8.3143 J/gmol/K), 
species molecular weight and mass fraction : 

n 

r> con 

~ (11.13) 

And the average molecular weight W m is defined as: 

1 ™ y. 

( IL14 ) 

11*2.3 Turbulence Model 

The two-equation effective diffusivity model is used to represent the turbulent char- 
acteristics. In the eddy diffusivity models, the turbulent fluxes, ^ and ^ , are 
related to the mean flow gradients through the assumption of an isotropic eddy 
viscosity and a constant turbulent Prandtl or Schmidt number: 


pu'iu'i 


f dUi 

- 


dUj , 2 

^-)+3 SiApk + vt 


dU k 

dx k ] 


( 11 . 15 ) 


pu'ift = 


Pt d^_ 
<j t dxi 


( 11 . 16 ) 


The eddy viscosity (/*,) appearing in equations (11.15) and (11.16) is defined in 
terms of a characteristic turbulence length scale(£ 3 / 2 / e ) and a velocity scale (jfeV*), 
so that p t is given by 


Pt - C^p— ( 11 . 17 ) 

The turbulent kinetic energy, k, and its dissipation rate, e, can be modeled from 
the turbulent transport equations: 
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dpk 

~dt 


d 


A »TJ.l-\ 


d 


IU 


file 


( 11 . 18 ) 


d Ji + JL iP u j( ) = /-(m + -)! 1 + - c,pt) + s,j ( 11 . 19 ) 

dt dxj KH 3 J tr/dxj k 

where the production term ptG takes the form 


PtG 


-rjdUi im dp dP 

pu ' u i dxj p 2 dx : dx } 

dU x dUj. ( dUi ^ dp dP_ 
dxj + dxi dxj p 2 dxj dxj 


- + /*< 


dlA.dUj 

dx r dx-j 


( 11 . 20 ) 


The last term of equation (11.20) will be zero for incompressible flow due to continuity 


equation. 


Table II. 1: Turbulence model constants 


C* = 1.45 Ci = 1.92 C„ = 0.09 a k = 1.0 ct £ = 1.3 a t = 0.9 


These are k — t equations [70] with some added particle modulation terms. Here, 
term involving in equations (11.20) is inserted to account for variable-density 
effects [57]. These terms originally come from the pressure- velocity correlation in the 
Reynolds stress equation. For reacting flows, these terms should account partially 
for the expansion effect on the flow field due to heat release from combustion. 

The added terms S k ,i and S c ,i accounting for the influence of particles on the 
turbulence structure will be discussed later. 


II.2.4 Combustion Model 


i 


Combustion of liquid fuel droplets in a spray is governed by the diffusion of fuel 
vapor and oxidizer species. Both premixed and diffusion flame theories can be used 
for spray combustion processes. Sometimes premixed-flame theories can be applied 
to some cases, where very small droplets of a high- volatility fuel may be completely 
evaporated in the heat-up processes. 

On the other hand, diffusion flame theories can be applied to many practical spray 
combustion processes. For these cases the fuel vapor evaporated from the droplet 
surface has to mix with the ambient oxidizer before chemical reaction can occur. In 
this study, it is assumed that liquid fuel droplets act as distributed sources of fuel 
which evaporate to form a cloud of vapor and the combustion process in spray flames 
can be treated as turbulent gaseous diffusion flames. Experimental evidence for this 
assumption can be found in [117]. An idealized approach for physically-controlled 
diffusion flames is to invoke a fast- chemistry assumption which the chemistry is 
sufficiently fast and intermediate species do not play a significant role. In the turbu- 
lent diffusion flame model, the influence of turbulence on combustion is taken into 
account by relating the fluctuations of mass fractions. This implies that fuel and 
oxidizer can coexist in the same place but at a different time. 

A modified eddy breakup model [76] is incorporated in the present study. Using 
this model, the reaction rate is determined as follows: in an irreversible single-step 
chemical reaction, the mixing-controlled reaction rate [76] is given by 


Rmix — A m i x p—rnin(Yf , — — ) 
k s 


( 11 . 21 ) 


where A miI - 4 is a model constant; s is the stoichiometric oxidant/fuel ratio; Y, 
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and Y 02 are the mass fractions of the fuel and the oxidizer. To account for the 
ignition delay time, the chemical kinetics need to be considered. The chemically 
controlled reaction rate, R c / ie , is given by the usual Arrhenius formula [119]. 

**• = O'- 22 ) 

The reaction rate, Rj u is determined from either of the mixing rate of the reactants 
or the chemical reaction rate, whichever slower. 


Rfu — TTllTli^Rmix i Rche'} 


(11.23) 


The irreversible single-step reaction of the hydrocarbon-air mixtures is expressed as 
follows: 


C x H y + (x + |)(0 2 + nN 2 ) -*• xC0 2 + V -H 2 0 + (x + |)nJV 2 (11.24) 

Here, n is 3.76 for air. In the given reaction process, five species ( fuel, 0 2 , N 2 , C0 2 , 
and H 2 0) are participating the mixture composition. Once the mass fraction of 
fuel and oxidizer have been determined from the solutions of the transport equa- 
tions, the mass fraction of the remaining species can be obtained from the following 
stoichiometric relations. 


Yh 2 o 
Yco 2 
V'n 2 
where 
K i 

I<2 


I< 2 ( 1 - K\Yo 2 - Yfu) 

I< 3 Yh 2 o 

1 — ( Yh 2 o + Yco 2 + Yq 2 + Yj u ) 


1 + n 


W N2 

Wo 2 


\W H2 q 

\Wh 2 o + (x + \)nW n2 + xW C o 2 


(11.25) 

(11.26) 
(11.27) 


(11.28) 


(11.29) 
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K 3 


xW C p 2 

2^f hO 


(11.30) 


II.3 Basic Lagrangian Equations 


In this study, the spray is described by a discrete particle method formulated on a 
Lagrangian frame. This is essentially a statistical approach and requires tracking 
a sufficiently large number of computational particles. Each computational particle 
represents a number of droplets having equal location, velocity, size, and tempera- 
ture. The particle characteristics are governed by equation of motion, equation of 
evaporation and equation of heat transfer. These ordinary differential equations will 
be integrated along the particle trajectories. 

II. 3.1 Droplet Motion Equation 


The equation of motion of a particle within a fluid continuum was originally derived 
by Basset, Boussinesq and Oseen for a fluid at rest, hence the B-B-0 equation, and 
extended by Tchen [114] to the case of a fluid moving with variable velocity. 


* , 3 

Pp ~dt 


\<PC D \V, - u,|(o, - u.) - ^ - £) 

3 t t dz^ 

+ fVftA (11.31) 


where t 0 is the particle starting time; the sub index p refers to the particle; ... and 
Ui are the velocities of the particle and the surrounding fluid respectively. The fluid 
velocity m should be defined at a distance far enough to the discrete particle not to 
be disturbed by the relative motion of the particle. 
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The term on the left hand side of equation(II.31) is the inertial force of the 
particle. The first term on right hand side of equation(II.31) is the drag force on 
the particle. The second term is the force on the particle due to the static pressure 
gradient in the flow field. The third term is the force due to the inertia of fluid 
displaced by the particle motion, or virtual mass term. The fourth term is the 
so-called ” Basset” term, which comes from the effect of the deviation of the flow 
pattern from steady state. The last term F bi represents the body force terms such 

as the gravity force and the centrifugal force. 

For the particular case of the motion of solid or liquid spheres in gas phase, where 
the particle density is much higher than the gas density, effects of static pressure 
gradient, virtual mass term and Basset force can be neglected. The drag term is 
treated empirically, assuming quasi-steady flow for single spherical particle. 

The approximate form of the B-B-0 equation of motion is 

dv i _ Uj ± ~ ^ + F bl (11.32) 

dt t 

In equation (11.32), F bl represents the body force terms acting on the particle per 
unit mass in i-direction, and u\ is included to represent gas fluctuating velocity effect. 
The particle relaxation time r can be expressed as: 

t -i = IL Cd |C/i +«/-*! (11.33) 

8 r p 


Cd is the drag coefficient given by 


Cd = 


f g-(l + \ReJ) if Re„ < 1000 


0.424 


if Re p > 1000 


(11.34) 


In which the particle Reynolds number Re p is defined as 

| Vi + Uj — Vjlpdy 


Re p = 


(11.35) 
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where p is the gas phase laminar viscosity. 

The particle position for each group can be obtained by integrating: 

dx{ 

IT = Vl 

II. 3. 2 Droplet Evaporation Model 


(11.36) 


Droplet size and temperature are governed by the mass and energy conservation law 
for each droplet. They are: 

d-T'p _ Th ey 

dt ^r^pd (11.37) 

and 

dT d _ Q l 

dt m p C p4 ( IU8 ) 

In equation (11.37), the droplet evaporation rate is given by the Frossling corre- 
lation [37] : 

m ev = 2ird p (pD)(l -f 0.3J?e p 2Sc p 3)/n(l _g m ) (11.39) 

In equation (2.16), the droplet temperature T d , which is assumed to be constant 
within the droplet, is found by using the heat energy Q L : 


Ql = 47rr p 2 Q c — m ev L 


(11.40) 


where L is the latent heat of vaporization, and Q c is the heat conduction rate to the 
droplet surface per unit area given by the Ranz-Marshall correlation [33] 


Qc = ^ (1 + 0.3Re p ^Pr p l)Ml±^l 




(11.41) 


The Schmidt number 5c p , Prandtl number Pr p , and mass transfer number B„ 
defined respectively as 


are 


Sc p 

Pr 

r 1 p 


JL 

P D 

p.C p 


I< 


(11.42) 

(11.43) 



and 
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y; - Voo v _ pj_ 
- t^y 7' y °° ~ P 


(11.44) 


The values of thermodynamical properties of gas such as K, C v , D etc. are 
highly dependent on the temperature and fuel vapor mass fraction at which they are 
evaluated. A ” one-third rule” [32] that utilizes a reference temperature equal to the 
droplet surface temperature plus one-third of the difference between the surrounding 
gas and droplet surface temperature is used. The same procedure is applied to the 
reference value for the fuel vapor mass fraction, in which Y s is obtained from 

= <"•«> 


Here Y s and P v are the mass fraction and the fuel vapor pressure at the droplet 
surface, and W } and W m are the molecular weights of fuel and mixture, respectively. 
For a given T d , P v is estimated from the J AN AF data bank [116]. 

The above Frossling correlation is valid only for dilute spray evaporation process. 
Radiative heat transfer and near critical and supercritical behavior are not included 
in this model. Also, more sophisticated evaporation models, such as Bellan s group 
evaporation model [10, 11, 13, 12, 14] or Chiu’s group combustion model [26, 53, 52, 
51, 25] should be incorporated in future studies to consider the dense spray effect. 


II. 3. 3 Droplet Size Distribution Model 


In a spray combustion chamber, sprays are generated by atomizers. These practical 
sprays generally consist of a series of non-uniform size droplets, or a spectrum of 
droplet sizes. Such a distribution of drop size varies under different liquid injections 
and operating conditions. Many experimental studies [72] were carried out to provide 
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Figure II. 1: Drop size distribution curves for a water spray in an air stream 


correlations for engineering designers and researches. A lot of mathematical models 
can also be used without the detail distribution information. 


Overall drop size characteristics are represented by their distribution curves 
which related to the cumulative percentage of droplet number (N), surface area 


(A), or volume (Q) as a function of droplet size (D), as shown in Figure II.l [79]. 

In the analysis of many spray dynamics, evaporating or combustion processes, 
only an average droplet size is desirable instead of the complete droplet size distri- 
bution. This droplet size is usually chosen as a mean or median diameter. General 
definition of a mean droplet diameter D ab raised to a power (a-b) is 


/ n _ So D a f(D)dD 
ab) T~DW)dD 


( 11 . 46 ) 


where 


f(D ^75 


( 11 . 47 ) 
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is a number distribution function, dN and N are the number of droplets in the size 
range from D to D+dD and the total droplet number respectively. Its corresponding 
cumulative function F(D) is defined as 

F{D) = f f(D)dD (11.48) 

Jo 

The commonly used mean diameter for approximate analysis of spray evaporation 
and combustion as an equivalent monodisperse spray is the so-call Sauter Mean 
Diameter” (SMD or D 32 ), which is defined as the ratio of entire spray’s volume to 
its surface area. 

The dilute spray model assumes that the fuel is injected into the combustion 
chamber as a fully atomized spray which consists of spherical droplets. The droplet- 
size distribution within the spray is represented by a finite number of size ranges. 
Several empirical distribution functions have been proposed to characterize the dis- 
tribution of drop sizes in a spray [72]. None of these is universally better than others. 
Model constants are adjusted to match the given set of data. 

X -Squared Distribution 

The normalized number distribution function for x -squared distribution is given: 

f{R) = -=4 exp(-=) l 11 - 49 ) 

6 R R 

where R is the number- averaged drop radius, which for this distribution is related 
to Sauter Mean Radius (SMR) by 

R = ]-SMR (11-50) 

o 

The corresponding cumulative distribution function is 

, i? xr R 1,R, 2 1/7? g. 

F(D) = l-exrf-jHl + j+jty'+jOlj)] 


(11.51) 
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This correlation needs only one constant (SMD or SMR) and is used in KIVA-II 
code [4], 

Nukiyama-Tanasawa Distribution 

f - (h.52) 

where a, 0, A, and B are experimentally determined constants. This expression is 
relative simple and adequately describes the actual distribution. 

Rosin-Rammler Distribution 

1 - Q = exp[-{j)<>] (H.53) 

or differential expression: 

dQ qD g ~ l D 

15 = ~xr ex P^x )q] ( IL54 ) 

x rn i, 

SMD ~ r(1 “ (n.55) 

where Q is the fraction of the total volume contained in drops of diameter less than 
D, and X and q are correlation constants. The relationship between X, SMD and 
q is established through the Gamma function T [72]. At present, this is the most 
widely used expression for spray drop size distribution. 

II. 3. 4 Particle Turbulent Dispersion 

Numerical modeling of particle turbulent dispersion based on the Lagrangian parti- 
cle tracking method was first proposed by Dukowicz [31] using a stochastic method, 
where particle turbulence was modeled by arbitrarily assuming gas turbulent kinetic 
energy and particle-eddy interaction time. It was further developed by Gosman and 
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Ioannides [39] and Shuen et al. [109] to include k - e turbulence model and estimate 
gas turbulent kinetic energy and eddy life-time. Main differences in the implemen- 
tations are the methods used to specify turbulence eddy properties and the methods 
for choosing the interaction time of a particle with a particular eddy. A two equation 
turbulence model, k - e, is used to characterize flow field turbulence quantities, such 
as turbulence fluctuation, eddy life time and length scale. Turbulence effects on 
particles are modeled by adding to the mean gas velocity U t a fluctuating velocity u\ 
when tracking particles through a continuous succession of turbulent eddies. Theo- 
retically, this simulation requires knowledge of the full time history of the turbulent 
flow by direct solution of unaveraged Navier-Stokes equation. Since this is impos- 
sible at present for most of flows of practical interest, the turbulence is simulated 
by means of a stochastic process or the Monte Carlo method. The instantaneous 
velocities for the gas phase are given by {/, + u\. Assuming an isotropic turbulence, 
each component of u\ is randomly chosen from a Gaussian distribution with standard 
deviation yj\k, where k is the specific turbulent kinetic energy of the gas phase in 
the computational cell in which particle is located. This assumption represents a 
significant simplification for actual complex turbulent flows where anisotropic tur- 
bulence will be expected, thus requiring a non-Gaussian distribution of turbulence 
intensity in three dimensional space domain. 

Adding fluctuating velocity u' to droplet motion equation (11.32), we have 

_ Ui + ~ + Fb t (H-56) 

dt t 

This equation represents particle-eddy interaction along its trajectory. The sum 
Ui+Ui is the gas phase velocity that transfers momentum to the particle which ’’sees” 
the eddy. It is this fluctuating velocity u/ produces particle turbulent dispersion, 
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Which keeps a piecewise constant function of time, changing discontinuous^ after 
passage of one particle-eddy interaction. The interaction time depends on the eddy 
life time (, and the particle transit time f „ within each eddy. 

To derive the eddy life time t„ we can choose a particle small enough that will 

fluctuate following the gas fluctuating velocity < We can expect in this particular 
case, 


= (11.57) 

where v' is the particle fluctuating velocity. Recall the assumption of the gas fluctu- 
ating velocity u,-, also v[ this time, obeying a Gaussian distribution with a standard 
deviation a u > 



After one particle-eddy interaction the deviation of the dispersed particle position 
will be 


Gx> = CT u 't e 

A distribution that follows a diffusion law has the following relation, 


(11.59) 


o x , - (11.60) 

In the limit of small particles, we expect that the particles will fully follow the 
fluid motion and their diffusivity D should be same as the eddy kinetic viscosity 
which means gas-phase momentum diffusion, 

D = v t = cjj. ( 11 . 61 ) 

The last relation coming from the k - e turbulence model used here. After some 
algebraic arrangements, we can obtain the eddy life time 

t e = 3 


(11.62) 
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Figure II.2: Shetch of the paxticle-eddy interaction 
Based on k — t turbulence model, an eddy dissipation length scale is proportional 
to fct jt and it has been used as the following 

l e = 1.65C£^- ( n - 63 ) 


where the proportional constant is chosen by fitting the experimental data of Snyder 
and Lumley [113]. 

In this study, the particle-eddy interaction time U n t is controlled by the eddy 
life time U and length scale l t . If a particle interacts with the eddy over a time t e 
or the displacement from the moving eddy center is larger than l e , the particle will 
interact with a new eddy and a new gas fluctuating velocity u\ will be generated as 
mentioned before. This process is shown in Figure II.2. 

Here, we do not estimate t int as suggested by Shuen et al. [31] in which t tT was 
calculated from a linearized equation (11.56) by keeping r a constant and without 
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the body force effect. Instead, we follow a similar stochastic procedure suggested by 
Nichols [86] and trace particle trajectories as time progresses. The turbulent fluc- 
tuating velocities u\ at three directions are generated with a Gaussian distribution 
when an eddy is created and remain constant in the particle-eddy interaction period. 
At each time step, the summations of the interaction time and the distance between 
eddy center and particle are stored at two arraies and compared with the eddy life 
time t e and length scale l e , which are calculated from equations (11.62, 11.63) based 
on the local gas-phase turbulence properties. If either of the eddy life time and 
length scale is exceeded, a new eddy will be generated at the particle’s location and 
the particle begins to interact with this new eddy. 

Figure II. 2 graphically shows an example of the particle-eddy interaction process. 
The particle first interacts with eddy e 2 and begins moving with the eddy center. The 
eddy time scale and length scale are t e j and l e i respectively. With time progressing, 
the particle remains in this eddy until the interaction time A tj exceeds t e i regardless 
^ 1,2 < l ei at that time. The particle begins interaction with a new eddy e 2 . After two 
interactions, the distance / 2 , 2 between the particle and eddy e 2 exceeds the length 
scale / e2 of this eddy and a new eddy e 3 is generated again. 

This method has the flexibility of taking into account both the gravity effect 
(crossing trajectory effect) and the non-Stokesian drag law and gives more satisfac- 
tory results for medium and heavy particle dispersions comparing with the experi- 
mental data of Snyder and Lumley [113]. 

The above procedure for solving equation(II.56) requires the calculation time step 
At being smaller than the particle-eddy interaction time otherwise the particle 
will ’’see” more than one turbulent fluctuating velocity u\ on the current cycle. 
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Possibly, we can use a smaller time step or subcycle time step for the particle phase. 
These methods are not computationally efficient, however, when t int is much smaller 
than At. O’Rourke [88] derived the particle velocity and position changes based 
on the assumptions of the Gaussian distribution of turbulent fluctuating velocity 
it' and constant particle relaxation time r, turbulent kinetic energy k and particle- 
eddy interaction time t int for the given particle on the current time step. O’Rourke’s 
approach [88] is inaccurate for the effects of the fluctuating velocity u\ on heat and 
mass exchange and droplet breakup and collision processes. 

The assumption of a linear drag law simplifies the analysis, therefore each com- 
ponent of the velocity and position changes can be treated independently. It has 
been shown[88] that for each component the distributions are Gaussian when the gas 
phase turbulent fluctuating velocity is Gaussian. The particle velocity distribution 
has variance 

o\, = j e3 T( ^ T ;j [l - exp{-2At/r)]al, (11.64) 

1 + exp(-t int /T) 
and its position distribution has variance 

0.2 

a\, = {tint At - 2t int T[l - exp(-tint/T)] + -yT 2 }a 2 , (11.65) 

a u' 

where <r 2 , = | k. 

When St > t tnt the particle velocity and position are updated using 


v’-v? ( Ui-v *) Sv'i 

~ET = + Fb - + Tt 


( 11 . 66 ) 


and 




( 11 . 67 ) 


where Sx' and Sv' are the particle turbulent position and velocity changes. 6x' is 
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calculated from 

bx = t per Sv' + Sx' b (11.68) 

Each component of 8v' and 8x' b is chosen from a Gaussian distribution with variance 
<r v , and <r x , — t per cr v , respectively. t per is a so-called turbulence persistence time and 
enters due to the independent distribution of turbulent particle velocity and position 
changes. 

*in<[l ~ exp(-At/r)] - ^pr~r 

W — ~ — (11.69) 

t 

II. 3. 5 PDF of Particle Turbulent Dispersion 


Assuming isotropic turbulence properties and a Gaussian distribution for gas turbu- 
lent kinetic energy, particle fluctuating velocity and position also obey this Gaussian 
distribution as seen from the above discussion. The particle PDF (Probability Den- 
sity Function) in two dimension form is 


p(*»y) 



2 a 2 ' 


(11.70) 


or in cylindrical coordinate 


pM) = 2 ^ eIp( -£ 


(11.71) 


where particles disperse from a point source (0,0). 

When calculations are performed in a cylindrical coordinate with constant prop- 
erties in tangential direction, any particle dispersion from a point (x 0 , y 0 ) in a trans- 
verse plane represent one from a circle with radius r 0 , where r 0 = yAo + Vo- This 
PDF can be obtained by integration along the circle (see Appendix A), 
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This PDF p(r, 9) represents particle probability in an area rd6 dr and is independent 
of 6. If we integrate the above equation over 9 domain, we will obtain particle 
probability in rdr , 


,, exp[-(r 2 + r5)/(2<r J )]^©)" 
PM = ^ ^ n\ 

n=0 


(11.73) 


In stochastic modeling of particle dispersion in two dimensional flows, usually we 
track particles in both x- and y-directions. 

dx)x x "h P 


— — = h F bx 

dt t 

dt T 

v z = 0 


(11.74) 

(11.75) 

(11.76) 


The equations (11.74, 11.75, 11.76) are true for two dimensional plane flows and incor- 
rect for axisymmetric flows due to the neglect of turbulent fluctuation in z-direction. 
To account for this three dimensional fluctuating phenomena, we solve particle mo- 


mentum equation in z-direction instead of keeping it as zero for axisymmetric flows, 


dv z u' z — v z 
dt t 


(11.77) 


and particle trajectory 



(11.78) 


Note that mean gas velocity and body force have been set to zero. Particle 
position in radial direction is calculated as 


r — 


(11.79) 


To test the present stochastic procedure, 10,000 particles are calculated at three 
different locations for each different circle, as shown in Figure II. 3, where particles 
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xi.o: r article dispersion 


» Vil 


begin to disperse. Particle distribution variance is assumed to be 1, which can be 
considered as particle position being normalized by the variance. Four circle radii 
(0.5, 1, 2 and 3) have been used to cover a large variation of radius-variance ratio. 

We can see from Figures II.4 and II.5 that this stochastic procedure results 
an equivalent particle distribution as described by equation (11.73) for an isotropic 
particle turbulent dispersion, and the distribution function is independent of the 
starting locations. Figures II.4 and II.5 also show the simplified PDF of Litchford 
and Jeng [75], that will be discussed in the next section, is in reasonable agreement 

wtth the present exact solution and good agreement as expected is reached for small 
and large radius-variance ratios. 



Figure II.4: PDF of particle dispersion from a circle (a) 







Figure H.5: PDF of particle dispersion from a circle (b) 
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II. 3. 6 Parcel PDF Model 


In this study, the spray is described by a discrete particle method formulated on 
a Lagrangian frame. The turbulence effects on droplet dispersion are simulated 
by a Monte Carlo method in the sense that a fluctuat.ng velocrty u' k , where each 
component of u\ is randomly chosen from a Gaussian distribution with standard 
deviation ^p, is added to the mean gas velocity. Thus the turbulence is assumed 
to be isotropic. This type of simulation for the turbulent dispersion of droplets has 
been extensively used previously [39, 36, 110] for statistically stationary turbulent 
dispersed flows and is called stochastic separated flow (SSF) model. In previous 
SSF model, computational particle position is characterized by a delta function in 
space. This computational particle represents a group of real physical particles 
which are assumed to stay at the center of the cluster. It does not consider the 
particle turbulent dispersion or probability function distribution within the cluster. 
When an insufficient number of particles are tracked and the two phase interaction 
source terms are evaluated, shot noise or over predicted flow field fluctuations will 
be generated in computational space and time domain. The shot noise will also 
increase in fine grid calculation. One of the important effects may occur in the spray 


combustion instability analysis when feedback responses to local transient flow field 
perturbation are over disturbed. Therefore, statistically meaningful calculations 
will require a large number of particle trajectories and, consequently, very long 
computational time or reduced spatial resolution of the analysis. 

To account for droplet turbulent dispersion, we follow the concept of Litchford 
and Jeng [75], and Zhou and Yao [124] of combining a normal (Gaussian) probability 
distribution for each computational particle. The instantaneous location of each 
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computational particle is calculated by a stochastic Lagrangian tracking scheme. 
The governing equation for each computational particle is 


dvk 

dt 


Uk - Vk 
Tk 


+ F b k 


dxk 



( 11 . 80 ) 

(11.81) 


with 


T, 1 = 


3 P Cd | . 

u r V ut ~ M 


(11.82) 


When u and v are taken as the instantaneous properties, the location calculated 
by the above equations only represents the mean of each particle’s corresponding 
probability function. The variance of each parcel PDF has to be calculated and the 


combined PDFs then represent the statistical distribution of particles with turbulent 


dispersion effects. To estimate the variance of the parcel PDF due to the turbulent 
particle dispersion, the turbulence-induced displacement and velocity can be splitted 
from equations (11.81, 11.82): 


dv\ _ tdfc - v' k 

dt Tk 

dx' k , 


(11.83) 

(11.84) 


With the isotropic turbulence assumption, each component of u'k is randomly 
chosen from a Gaussian distribution with standard deviation u'k rms = \J\k- We 
first choose A tk t as the time step of the particle i th interaction within the k tfl eddy, 
which is smaller than its eddy life time t e k , and integrate equations (11.83, 11.84) in 
that sequence to update particle fluctuating locations and velocities. 


x 'ki ~ u'krmsAtki + 1) - V, 1 k rms j) ( 1 - ex p(- 


Af fc; 


ki 


A-b-i) 


)) 


(11.85) 
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v'ki = u'krms + ( u ^(t-l) — u ' krms) ex P{ ~ ) (11.86) 

T *(*'- 1) 

We then sum up the m steps for which the particle fully interact with the k th 
eddy, 

m 

At** = At* = t e k (11.87) 

i=i 

The change of variance of a computational particle PDF within the k th eddy is 
represented by a characteristic mean squared dispersion in the form: 

m 2 

^ = ^-,’ + £ 4 ) ( 11 . 88 ) 

i=l 

In equation(II.88), (Tk-i is the existing variance of the particle PDF at the beginning 
of the interaction within the k th eddy. Since the time step within each turbulent eddy 
is fixed, the number of interaction within the eddy, m, varies across the calculation 
domain, the choice of time step At*,, and the related issues are discussed in detail 
in [74]. Figure II. 6 describes this eddy interaction with the particles. 

The present procedure is easy to program and requires less computer memory. 
For each computational particle, we just need to store x'ki , u'krmsi v 'ki, and cq, 2 . 
This procedure when implemented in the current time-marching numerical method 
is somewhat different from the method of Litchford and Jeng [75] in which the calcu- 
lation of the current variance of each particle PDF is summed over the entire history 
of the effective time constants. In their recent study, truncation of unnecessary time 
history terms and the associated errors were discussed and additional computational 
efficiency was obtained [74]. 

When convoluting PDF for a group of computational particles with their trajec- 
tories calculated by SSF model, the variances of equation (11.88) must be normalized 
according to the total number of particles. The normalized particle variance can be 
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Figure II. 6: Eddy interaction with the particles 


written as 


a yk = K-^L 


Her,, ^ represents the statistical uncertainty in the mean particle position, K is 
the correction factor to account for undersaving, and * is the total number of 
computational particles. When symmetry and motive boundary condition exist in 
the calculation domain, a cumulative PDF distribution at any point in coordinate 

y, which is the distance from the particle to the axis or the reflective boundary, can 
be defined for plane and axial svmmptnr . . . 


Pamcle to the «*■ or the reflective boundary, can 
symmetric coordinates respectively. 


Plane Symmetric 


• P(y) = r ^Jr~ 

J-y V2ircr yk 

Here, y p is the instantaneous location of cc 


expt-^dy 

yk 


(11.90) 


computational particles. After integration, 
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the symmetric cumulative distribution function takes the form 

,y -y P s . r/V + yp 


P(y) = 0.5[er/(- 7 =7-^) + er/(^ i 

where the error function is defined as 


&yk 


)] 


9 f x 

erf(x) = -7= / exp(-C 2 )^C 


The corresponding PDF is 


0.5 


(y — Vp )2 i . el p|_ 

~ 2 » ir i+e vl 2 <7 




(11.91) 


(11.92) 


(11.93) 


Axial Symmetric 

Particle PDF distribution at this case should be equation (11.73), but it is too com- 
plex to integrate. A simplified axisymmetric cumulative distribution function intro- 

duced by Litchford and Jeng [75] takes the form 

F{r) (11.94) 


P(r) = 


F(r — ► oo) 


where 


— . rl , (r-r P ) 2 , rilllz)!]} 

F{r) = V2^a r {2exp(-^)-exp[ J exp[ ^ JJ 

r f( r ~ rp ) er f( - ~^~ — ) + 2er/(— p — )] ( IL95 ) 

+ 7rrp[er/( 7!d7 ) n V2a r j V2d/ 


and 


F(r -> oo) = 2v / 27rd r exp (- ■—) + 27rr p er/( ^ 


-) 

<7r 


(11.96) 


In accordance with the approach of Litchford and Jeng [75], when the mean 
positions of computational particles are calculated by the deterministic tracking (»* 
in equations (11.80 - 11.82) is the mean gas velocity), this approach is described 
as the deterministic dispersion width transport(DDWT) model. For tracking with 
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stochastic sampling of gas-phase turbulent velocity fluctuations (u k in equations 
(11.80 - 11.82) is the instantaneous gas velocity), the approach is described as a 
stochastic dispersion width transport (SDWT) model. 

II. 3. 7 Turbulence Modulation Model 

Turbulence modulation effect appears with the presence of the dispersed particle 
phase. This effect comes into the governing equations (11.18 and 11.19) of turbulent 
kinetic energy and its dissipation rate through the source terms S k ,i and S t /. Shuen 
[107] derived the expressions based on the gas-phase momentum equation (II. 2) with 
the interphase source term S u ij. These terms can be written as 


S k ,i = u t S u ij — U t S., 


? uij 


(11.97) 


5. 


e,l 


= 2 ii 


duj'dS^j 
dxj dxj 


(11.98) 


Within the framework of discrete particle stochastic approach using Lagrangian 
tracking, u' t follows the Gaussian distribution and the instantaneous properties of 
two-phase interaction force S ui j takes the form 


1 i v r t 

Sui ’‘ = dV ~ ™pN P ( - t + Ul ~ Vt ) p ] (11.99) 

p = i 7 

where u,- is the instantaneous particle velocity, r is the particle relaxation time, m v is 
the droplet mass and m ev is the droplet evaporation rate as defined before. Equation 
(11.97) can also be written as 


S k ,i = u'iSuij (II. 100) 

Thus can be calculated directly without modeling in equation (11.18). This is 

the approach used in [110, 4]. For e equation, the modulation term must be modeled. 



32 


The proposed model of Amsden and O’Rourke [4] will be used here in association 
with the k equation. This model is assigned Model 2 in this study in which the 
extra term in e equation is 

St, I = c 3 ^s k ,, ( 11 . 101 ) 

with C 3 = 1.5. 

Substitution of equation (11.99) into (II. 100) and taking average, we obtain the 
modulation term 

1 » p — 

Su = 1U - n. p JV p u'(-i— -i) ] (11.102) 

dv * — ' t p 

p= 1 

Mostafa and Mongia [80] as well as Fashola and Chen [36] have proposed a 
simplified approach in which the interaction term S ut ,i is linearized followed by mul- 
tiplication of the fluctuation velocity u\. The turbulence modulation term then 
only involves the gas/droplet velocity fluctuation correlation u'(u' — t>') in equation 
(11.102) for non-evaporating sprays. This correlation is then modeled through gas 
kinetic energy k, and eddy and particle time scales, ti and t d , 

- u') = 2k f(U,t d ) (11.103) 

In Mostafa and Mongia’s model [80] ( Model 1 in this study), the correction 
function / takes the form 

/((„« = 1 (11.104) 

where ti is the gas phase Lagrangian time scale given by 

U = 0.35 - (11.105) 

and t d is same as the particle relaxation time r defined in equation (11.33). 
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Following approach of Chen and Wood [22] ( Model 3 i 
in their two-fluid model, we use the function 


— 1 — exp(— — ) 

U 


in this study) employed 


(11.106) 


where 


and 


t, = 0.5 - 

e 




The extra term in k equation 


18 p/j. 


(11.107) 


(11.108) 


spray, 


is summarized as followings for non-evaporating 


NP 3 

^4d p C ^ u ' ~ v, 'lp 2 Tr p 3N , ik f(ti,t d )/dV (11.109) 

The corresponding extra ternr in s elation is modeled as equation (1U01) 

with constant C 3 = 1.0. In Model 1 and 3, effects of different turbulence timescales 

wth respect to parttcle relaxation times are incorporated m the modulation terms 

Furthermore, this mode, simplifies elation of the dispersed-phase source terms in 
two-phase flows. 

Based on the above modulation models, the modulatton terms appearing in * 
equation, S w , have negat.ve values al, the time, hence the presence of panic, e-phase 
will damp the gas-phase kinetic energy and affect the turbulence structure. 

H‘3*8 Droplet Breakup Model 

The present study employs two breakup models including Reite's wave instability- 

model (95J and TAB (Taylor Analogy Breakup) model of O'Rourke and Amsden 
[89). 
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Reitz’s model 

Based on the concept that the atomization of the injected liquid and subsequent 
breakup of drops are indistinguishable processes within a dense spray N, «-U 
and Diwaker developed a droplet breakup mode,. This model was extended by Re.tz 
[ 95 ] with adding new parcels containing fine product drops to the computat.on as 
the blobs break up. Atomization is modeled using injected ’blobs-, which have the 
same sizes as the nozzle exit diameter. The breakup of the blobs and the resultmg 
drops are modeled using a wave stability analysis for liquid jets. The wavelength of 
fastest growing wave, A and the maximum wave growth rate, 0 can be determined 
by the curve-fitted formula which are obtain from the numerical solutions of the 

surface wave dispersion equation for a round jet. 

A (1 + 0.45Z o - 5 )(l + 0-40r°- 7 ) (11.110) 

Z~~I 9.02 ✓ - s\ n r T I , r . 1 fi7\0.6 


a 


(1 + 0.87 We*- 67 ) 0 - 6 


Q 




3 \ 0 5 0.34 + 0.38We 


1.5 

9 


(1 + Z){ 1 + 1.4T 0 - 6 ) 


where 

W = \U + u'-v\- gas-droplet relative velocity 

^ = liquid Ohnesorge number 

yf ei _ pjj°W 1 liquid Weber number 

y/ 6g= paW - — gas Weber number 

jjf — liquid Reynolds number 

T= Z WT g 

a - liquid jet radius ; a - droplet surface tension 


(II. HI) 
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The mean product drop size and drop breakup rate are given by 




Boh if B 0 A < a 

2iT7/ork\°- 3 3 


min 


{3na 2 w/2ny 
(3a 2 A/4) 0 - 33 


if B 0 A > a 


( 11 . 112 ) 


and 


da (a — r) 
dt r b 


(11.113) 


where 


r 6 = 3.726B\a/ ACl (11.114) 

Here, B 0 is 0.61 and B\ is the breakup time constant, a is the radius of the liquid 
jet or the blob. The secondary breakup is assumed to be governed by the same 
equations for the primary jet breakup. The finer drop parcel is generated when its 
mass reaches 20 % of the parent drop mass. The breakup constant Bj = 10 is used 
for atomization process, and as suggested by O’Rourke and Amsden [89] B x = 1.73 
is employed for droplet secondary breakup. 


TAB Model 

The TAB model of O’Rourke and Amsden [89] is based on an analogy between an 
oscillating and distorting droplet and a spring-mass system. The restoring force 
of the spring is analogous to the surface tension forces. The external force on the 
mass is analogous to the gas aerodynamic force. The damping forces due to liquid 
viscosity are introduced to this analogy. 

Based on the TAB model, the equation for the acceleration of the droplet dis- 
tortion parameter is 

= % P (U + u’ -v) 2 _ 8a(T d ) _5pi(T d )dy 
dt 2 3 pd r 2 p d r 3 ^ p d r 2 dt 


(11.115) 
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where the quantity y is proportional to the displacement of the droplet’s surface 
from its equilibrium position, divided by the droplet radius, a and pi are the liquid 
surface tension coefficient and viscosity respectively at the droplet temperature. 


After integrating equation (11.115), we have 

We t We 1 

y(t) = — + exp{--)[(y( 0) - — cos ut + — (y(0) + 


where 


We = 


p(U -f it' — v) 2 r 


^ — ) sin ujt] (11.116) 


(11.117) 


is the Weber number, 


is the viscous damping time, and 


_ 2 p d r 2 
- 
5 


u? = 8 


a 1 


(11.118) 


(11.119) 


is the square of the oscillation frequency. 

The droplet oscillation and breakup calculations require two normalized particle 
arrays(y, deformation and oscillation) which can be determined by the above 
equations. Droplet breakup occurs if and only if y(t ) is greater than unity. Occur- 
rence of droplet breakup, the Sauter mean radius(SMR), and oscillation velocity for 
the product drop depend on these two parameters and Weber number. The radius 
of the product drops is then chosen randomly from a x -squared distribution with 
calculated SMR, which is given by the following relation: 


SMR = 


1 _L 

3 ' 8 a V dt ) 


( 11 . 120 ) 


where v is the parent droplet radius and ot is the droplet suiface tension. Following 
breakup, the product drop has the same temperature with the parent drop, and its 
deformation and oscillating parameters are set to zero. 
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II.3.9 Droplet Collision Model 


The drop collision model suggested by O'Rourke (S7J is employed to calculate col- 
lision and coalescence among the dispersed liquid phase. The collision routine is 
operating for the pair of particles if, and only if, they are in the same computational 
cell. For the collision calculation, the drops associated with each computational 
parcel are considered to be uniformly drstributed throughout the computational cel. 

they are located. For all parcels in each computational cell, a collision fre- 
quency between the parcel (parcel.) of larger drop radius(r,) and the parcel (pa rce W 
of smaller drop rad.us(r 2 ) is obtained from the relationship in terms of the number 
of drops in parcel,, the relative velocity between parcel, and parcel 2 , the area based 
on r, + r 2 , and the volume of computational cell. Such a frequency a is expressed as 

U = dVZ 7r ( r l +r 2 ?\vi~ v 7\ ( 11 . 121 ) 

The probability P for » collisions is assumed to follow a Poisson distribution based 
on a collision frequency and the computational time step At. 


Pn = e-*2- 

n! 


where the mean value n = uAt and the collision frequency is defined 

( 11 . 121 ). 


( 11 . 122 ) 


m equation 


Using the probability informations, the collision impact parameters are stochas- 
tically calculated. If the collision impact parameter is less than a critical impact 
parameter, the outcome of every collision is coalescence. Otherwise, each collision 
is a grazing collision. The critical impact parameter depends on the drop radii, the 
relative velocity between drops, and the liquid surface tension coefficient. 
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Suppose the outcome of the collision is coalescence. For each collector drop, 
n droplets are subtracted from their associated parcel, and the size, velocity, and 
temperature of the collector parcel are appropriately modified. If there is not enough 
number of droplet to have n coalescences with each collector, then n is recalculated 
so that all Hi droplets coalesce, and the parcel associated with these droplets is 

destroyed. 

If the outcome of each collision is a grazing collision, only one collision is calcu- 
lated for each parcel. Grazing collisions usually occur between drops of nearly equal 
size and are calculated between N pairs of drops, where 

N = min(N^N?) ( IU23 ) 

The droplets maintain their S1 zes and temperatures but undergo velocity changes. 
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II. 4 Summary of Two-Phase Interaction Source 
Terms 


The two-phase interaction source terms in the governing equations can be sum- 


marized as : 



j NP 

S m J — ^ ^ ^ v Ta ev>v 

p=l 

(11.124) 

SuiJ — 

1 ^ Ui + u'i - v t . . 

— ^P[N p m ev , p (vi) p - m p N p { r )p] 

p=l 

(11.125) 

NP 

p=l 

1 o dTp dv{ . 

rn ev , p (hp — L + ) — m p N p (Cp t d ^ ^ v *)} 

(11.126) 

where 

du, _ + «<-”*■ +J p 6i 

(if p r 

(11.127) 

and the particle mass 

_ ^ 3 

TTlp j^^pPd 

(11.128) 


The droplet evaporation 


rate can be expressed as: 


4 ( r p +1 ) 3 — ( r p ) 3 


m, 


= Ti^pd- 


At 


(11.129) 


Turbulence modulation terms in k and e equations take the following forms, 


Model 2 (Amsden and O’Rourke [4] ): 


Sk,l — 


(11.130) 


i = 1.5 -S kt , 


(11.131) 


Model 1 (Mostafa and Mongia [80]) and Model 3 (Chen and Wood [22]): 

NP 


3 P_ 

i = 1 - ip 


Su = -Y it -C D \ui - Vi\, -7rr p 3 N p 2k f(t h t d )/dV 

’ ' 4 dr, 4 


(11.132) 
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Sc, I = l.o|s*,, 

The velocity fluctuation correlation function is defined in section (II. 3. 7). 

Here, dV denotes the volume of the computational cell and h p and 
droplet enthalpy and the latent heat of the droplet, respectively. 


(11.133) 


are the 


Ill 

Numerical Procedure 


Introducti 


ion 


The governing equations of the gas phase are solved using a contro,. volume based 
finite-difference method on an unsteady fashion. Spatial differences are formed on 
a curvilinear general coordinate with all gas field variables stored at the same grid 
point. Second order accurate central differencing scheme is used for the diffusion 
terms and high order Chakravarthy-Osher scheme [19] with damping is used for the 
convection terms. The transient solution is marched forward in a sequence of finite 
time increments. The implicitly coupled pressure and velocity equations are solved 
by the M - PISO algorithmic], with individual equations being solved by the 
conjugate gradient squared(CGS) [29] method. In the P,SO algorithm, each time 
step is divided into a one-predictor two-corrector sequence. 

The strong coupling terms between particle and gas are evaluated by the same 
time splitting technique. Implicit coupling procedures are used to treat momentum 
exchanges to avoid the the limitation of small time-steps. 

Accurate caicuiation of mass and heat transfer is achieved by automatic reduc- 
tions in the timestep when the exchange rate becomes large. For droplet/turbulence 
interaction calculations, integration time step is compared to the turbulent eddy life 
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time. If the time step is smaller than the eddy life time, a fluctuating component is 
added to the local mean gas velocity when calculating each particles mass, momen- 
tum, and energy exchange with the gas. If the time step exceeds the eddy life time, 
changes in droplet position and velocity due to turbulence are chosen randomly from 
the probability distributions for these changes as described by 0'Rourke[SS]. 

The unsteady solution procedure described above is different from the conven- 
tional PSIC (Particle Source In Cell) procedure [28] in which global iterations are 
required between two phases. The method used here is time-accurate and is very 
computational efficient. This unsteady procedure can also be used for steady state 
calculations where the statistically steady solutions are sought. Detailed descriptions 
of the current method are given in section (III.5). 

XII. 2 Statistical Particle Model 

In this model, spray is represented by discrete particles, rather than by continuous 
distributions. A finite number of computational particles are used to predict a sample 
of total population of particles. Using the statistical or Monte Carlo formulation, 
each computational particle is considered to represent a group of particles having 
the same characteristics such as number N„ size r p , velocity v p , location * p and 
temperature T d . The discrete particle distribution function / is used to approximate 

the continuous distribution, 

f = N ^ x ~ Xp ^ v ~ Vp ^ r ~ rp ^ T ~ Td ^ dV 

P=1 

Particle trajectories are integrated using its motion and momentum equations 
(11.32,11.36) and particles exchange mass, momentum and energy with the gas within 
the computational cell in which they are located. By re-arrangement of two phase 
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momentum interaction terms, this method consists of a fully-interacting combination 
of Eulerian fluid and Lagrangian particle calculations. The interaction calculations 
are performed simultaneously and eliminate global iteration for the two-phase mo- 
mentum exchange. This procedure improves numerical stability and is efficient for 
transient calculation. Dukowicz [31] firstly proposed this numerical procedure and 
used a random turbulence model to simulate unsteady spray jets. 

In the stochastic particle method, a Monte Carlo technique is employed for direct 
simulation of spray characteristics. Droplets are sampled randomly from assumed 
probability distribution functions that govern droplet properties and droplet be- 
havior subsequent to injection. The sampling procedures include injection droplet 
size and velocity, gas turbulent fluctuating velocity, droplet breakup and droplet 
collision. This procedure is summarized as following. 

In the direct simulation of Monte Carlo procedures, first we need to find a table 
of successive random fractions R s that are uniformly distributed between 0 and 1. 

This random number can be generated as a standard function on computers or from 
some algorithms. 

Suppose we are given the droplet number distribution function fir) correspond- 
ing to the random radius r (r, < r < r 2 ). The droplet number dN in the interval dr 
about the radius r will be 


dN = f(r) dr 


(1H.2) 


We define the random variable 


^ “ I n ^ ^ ' (III.3) 

and we know that dN = dr f . So, the number of droplets is uniformly distributed 
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with respect to the variable r/. If we normalize function r } by 

f 2 f(r')dr' (III.4) 

J T\ 

and define 

Rj= f f(r')dr'/ ! f{r')dr (HI.5) 

Jti J r\ 

Hence, R f will be distributed also from zero to one, and then we invert Rf to obtain 
r, which then will be distributed according to f(r). If the form of the distribution 
function f(r) is easy to integrate and its integration is easy to invert, we can perfoim 
this process analytically. Otherwise, we have to use a numerical method. 

The following examples illustrate how these processes work. Consider the droplet 
radius r distributed between 0 and a such that the probability of r is proportional 


to r. Now 


u. 

II 

(Hl-6) 

and 


r, = \\ir = \? 

(III.7) 

Normalized function Rf will be 


r 2 

Rf — ~2 

a 1 

(III.S) 

Its inversion is 


II 

S- 

(III.9) 

Unfortunately, the above procedure can only be used when it 

is possible to obtain 


an explicit function for r. 

Let us consider gas turbulent fluctuating velocity component u'. Due to the hy- 
pothesis of isotropic gas phase turbulence, u' is distributed according to the Gaussian 
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distribution 

f( u ) — ——=-exp( ) 

V2 xa y 2a 2 ’ 

where, a is the standard deviation given by: 

fT~ 

’ = V 3 A 

and k is the gas turbulent kinetic energy. Its integration will be 

R ‘ = L ^h cxp( ~^ du ' = 5 + rr/fuvysTJ) 

It follows that the error funct.on erf m us« be inverted to obtain u ' 


(III.10) 


(III.ll) 


(III. 12) 


u ' = V^kJZerf~\2R s - 1 ) 


(HI. 13) 


or 


- V4k/3s,gn(2Rj - \)erf-\\2R, _ Z |) (nU4) 

Where sign is a FORTRAN function which takes the sign of a variable. 

We store the values of the tnverse error function erf^R t) m an array ^ 
an intervals of 0.05 from R, = 0.0 to R, = 1.0. The maxima value of 
erf-'(l.O), is infinite and approximated to be er/-(0.99532) = 2.0 without sig- 
ffect on calculation of u , because values of uniform distribution random 
function Rj exceeding 0.99532 have very small probability. Values of erf^R,) at 
intermediate values of R, are obtained by linear interpolation. A new „ is only 
sampled once every particle-eddy interaction time t, M . 

HI.3 Splitting of Two-Phase Momentum Inter- 
action Source Term 


To improve the 
source term, S uit 


convergence and the numerical stability, the momentum interaction 
, can be treated implicitly. The particle momentum equation (11.32) 


46 


can be discretized as 


vj^^ = ui2±AzJ±l + Fbt 

At 


rTl + l 


(III. 15) 


and v t n+1 can be analytically solved from the above equation, 

; n + ([/n+1 + u *. + iq,jT n+1 )^r 


u” +1 = 


1 -1- 

i T* T n+1 


(III. 16) 


Thus, 


[ 7, n+1 + u' - 


Tl + 1 __ 


(V, 


■n+1 + u ')(j + -£L.) - [tff + (6/f +1 + Uj + ^ T " +1 ) T"tl 


1 4- -£A. 

r > T n+1 


Two-phase momentum o 
above relation (III. 17) 


U) 1+1 + u' ~ jwAt - v? (III. 17) 

“ 1 + ^ 

oupling term (11.125) can be obtained by substituting the 


, NP . , [/" +1 + u\ - u" +1 A 1 

sifj = 2v^ NpThev ' p ^ np ~ mpN ^~ ^ "p 

P=1 

= + E — (Vi) ' ■ TT^ (r ' 


dV 3 

P =i 


.(C / t n+1 +u' -FwAt-v")] 

NP 4 /vr (r n+1 ) 

_ j?^A{( r ") 3 (t, i ) w 

- av A - 1 p 


ciV Z-jr At 
p-i 


P 1 + 7 


3 

aT 

n+l 


At 


.(([/»+> + u ;-F l ,Ai)^ Tr + »”l) 

We then split into the following expression, 

S3 1 = -S pt U? +1 + Rr* 


(III. 18) 


(III. 19) 


Here, S p and R p are 


obtained by comparing equation(IlI.lS) and equation(HI.19), 


S p i - 


j_ " \,p d N r (O 

dV ^ 


P=1 


r n+1 1 + ^r 


(III. 20) 
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Rju = ~ i*pdN p 

+ (ni.21, 


P=1 


°-{(r n p ) 3 (v,) n - i5 +I ) 
1 + 


T^ + l 


T he parameters Sp and * are momen(um contro] vQ|ume quantjtie$ 

"* ° n aVailaWe Partide “ f “° - ~ ™* P- nt tnnestep, and * „ 
independent of directions. This splitting w a s firstly proposed by Dukowicz , for 

" octy to the two-phase conpling source ternr a „d enhances two-phase interactions 
n strong two-way coupling calculations, its efficiency and numerical stability have 
>een shown by Dukowics [31, and the present study. For transient calculations 

accuracy solutions. This splitting enhances diagonal dominance of the momentum 

IIL4 Dr ° plet Evaporation Calculation 

One of the advantages of hagrangian method for li q md phase ,s its feasibility to 
andle polydispersed spray. Evaporation rates of various sires of spray droplets^ 

Be calculated for each droplet, and then the twmphase mteractrons can be coupled bv 

source terms. Droplet radius and temperature changes are calculated by evaporating 

droplets sequentially at constant pressure Tn 

P e ' In hlghly eva P°ratmg processes, such the 

■gh ambient temperature or for more volatile fuel droplets will 

1 aropiets will evaporate very 

quickly and small time step will be reouired , 

q red to get reasonable numerical solution. 

To improve computational efficient n • 

efficiency, following the numerical treatment of K1VA-II 

W “ hybtid meth ° d iS ^oration calculation of droplet is 

implicit in its temperature hut explicit in the gas temperature and mass fraction. 
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Droplet evaporat.ons are emulated one by one along the.r trajectories when they 
travel through the god cells and gas phase temperature, pressure and mass tract, on 
at the grid point are used as ambient properties. This sequent, a, droplet evaporat.cn 
and explicitness can produce unphysical results and numerical instability 

an evaporation sub time step *. is estimated based on heat and mass transfer rates 
for each parcel of droplets and the evaporation calculation is subcycled with fl. to 

main time step At. 

The choice for «„ is based on the idea that the heat or mass transfer between a 

computational particle and its surrounding gas phase at the same grid point in one 

time step should not exceed some fraction of energy or mass available for transfer. 

This fraction is specified to be 0.5 here. This criterion can be expressed as 

pdV (III. 22) 

6tev ~ uSh 4t rr p iVp 

_ m 391 we have droplet size change equation 
Rearranging equation (il.oyj, 

dr 2 p (pD)Sh In(l + B m ) (III. 23) 


dt 


Pd 


and temperature change equation 

iTi K(T - - L (i>D)Sh ln(l + B "L> 

dt pd$ r p ( ~'pA 

number Nu and Sherwood number Sh are given by 


Where, Nusselt 


Nu = 2 + 0.6 Re p 2 Scd 3 


(III. 24) 


(111.25) 


and 


Sh = 2 + 0.6 Re p ^Pr d 3 


(III. 26) 
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First, droplet energy equation is solved implicitly to update its temperature. The 
finite-difference approximation to equation(III.24) is 


T*'-TZ K*{T-T?'){Nv) 


ylMl+B# 1 


f- 1 ~ L" (pD)" (S h)" ln(l + By") 


pdK r p) 2 cZd 


(III. 27) 


where the superscript v denotes the value of a quantity after v evaporation subcycles. 
The Spalding number or mass transfer number B m is calculated from the formula, 


B u . 

m i - t s (t; +1 ) 


(III. 28) 


Yoo in equation (III. 28) and T in equation (III. 27) are intermediate values of ambient 
vapor mass fraction and gas temperature of the grid cell in which the droplet is 
located. They were explicitly updated due to evaporation of droplets with subscripts 
less than current particle index p and evaporation of this droplet on subcycle less 
than v. The gas phase thermal conductivity K v and mass diffusivity ( pD) v are 
calculated using T and T%. The droplet Reynolds number Re p , Schmidt number 
Sep , and Prandtl number Pr p are calculated using r£, and the intermediate gas 
temperature T. The finite difference equation (III. 27) requires iterations since mass 
transfer number or droplet surface vapor mass fraction Y s is a nonlinear function of 
droplet temperature Tj. 

The droplet size change will be calculated following the implicit solution of equa- 
tion (III. 27) for T % +1 . The finite difference equation is 


K+'r - Kr {pdy „,Mi + K) + '"(i + ^ +1 ) 


(III. 29) 


If the newly calculated droplet size r£ +1 is smaller than zero or temperature T% +1 is 
greater than its critical temperature, the droplet radius will be set to zero and we 
will terminate the evaporation calculation for this droplet. 



The intermediate gas temperatures and vapor mass fractions are calculated ex- 
plicitly with droplet evaporation progressing. Three arrays related to gas mass, fuel 


vapor mass and gas enthalpy in grid cells are initialized as: 

Mij = pijdVij (III. 30) 

(M f h = (pf)ijdVij (III.31) 

(tfH)n = PijdVij h{j . (III.32) 

Then as each droplet evaporating the above arrays are updated as 

Mi, - Mi, - - {r ‘;,y\ N i (in-33) 

(Mf)i, <- ( M,)„ - - Kf]N r (III. 34) 

(MH)ij - (MH)i, - fh"*' - (r;fh;\N r (III. 35) 

V r p d [(rrf - (r;) 3 ]N r L-* 1 (III.36) 


for the droplet located in the cell (i,j). Finally, the intermediate vapor mass fraction 
and gas temperature are calculated, 


(n.)« 

(M,)„ 

Mi, j 

(III. 37) 


( MH)i, rrn 


T„ 

1 

1 

Pi-D 

§ 

+ 

Eh 

II 

(III. 38) 


These new intermediate values are then used for the next calculation of droplet 
temperature and radius changes. We can get droplet evaporation rates and their 
contributions to energy equation due to evaporation after evaporation calculation is 
completed. Droplet mass evaporated is 
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Energy contribution i 


is 


NP 


P =1 J F 

- EfV'f (’■’if - (C'f)N rL -» } (IIL40) 

In the above equation (III.40), £„ denotes that the suction is taken oven a, I of the 

■ub- time step Because the latent heat L- depends on the droplet temperature 

T d » we do not have such a relation: 

E rr<Wf - = t xpiKr;f _ (f .»«)3 ]A;i . +1 (m 4J) 

When impiemented in the two-phase calcuiation, the droplet mass evaporation 
term S mJ and energy exchange term S v are employed as: 


NP 


c 1 y — ' pdN v 


P=zj 


(III. 42) 


Sk,i = 


- E 


1 ly r 4 TK T 

-L r ^Ttpd^v 

dv - ( r ; +1 ) 3 (^)" +1 ] 

i^PdNt, 

a T~ f(r ^ 3 " K +1 ?]l u+ '} 


(III. 43) 


III.5 Two-Phase Numerical Model 

Due to the strong coupling and stiff source terms in the two-phase flow, successful 
numerical schemes must be stab.e, accurate and efficient- The spirits of Fssa’s [47J 

P,S0 alg ° rithm “ d DUW1C2 ’ S [31) -‘-‘-I and splitting techniques are com- 
bmed and a new two-phase coupling scheme is proposed in the following sections- 
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III. 5.1 The PISO Algorithm 

The PISO (Pressure-Implicit with Splitting of Operators) algorithm [47] is a non- 
iterative method for handling the coupling of the implicitly discrete time-dependent 
single phase fluid flow and heat transfer equations. The method is based on the use 
of pressure and velocity as primitive variables and hence can be applied for the so- 
lutions of both incompressible and compressible versions of the transport equatrons. 
The main feature of the technique is to split the solution process into a series of 
predictor and corrector steps whereby operations on pressure are decoupled from 
those on velocity at each time step. Due to the splitting, the set of equations can be 
solved sequentially. At each time-step, time accurate solutions can be obtained by 
prescribed predictor and corrector steps. The accuracy of this splitting procedure is 
based on a linearized form of the discretized equations, and the analysis indicates 
that the numerical solution differs from the exact solution of the difference equa- 
tions by terms proportional to the powers of the time-step- size. By virtue of this, 
it is possible to dispense with iterations and work faster for transient flow calcula- 
tions. This efficient implicit scheme retains simplicity of implementation relative to 
block simultaneous methods which require much more computer memory and are 
not flexible to include extra physical models. 

III. 5. 2 M-PISO and Gas Phase Solver 

As pointed out by Jiang [55], the original PISO algorithm suffers some problems in 
splitting procedure and is not applicable to supersonic flows. A newly modified M - 
PISO algorithm proposed by Jiang et al. [55, 56] has been successfully earned out 
for incompressible, transonic and supersonic flow calculations. The benchmark test 


53 


cases[55] include incompressible channel flow, driven cavity flow, vortex shedding, 
both laminar and turbulent backward facing step flows, and compressible flows over a 
bump and a spherical cylinder, nozzle flow and multiple shock reflections. M-PISO 
algorithm is employed and improved in the present gas phase solver. The numerical 
details can be found in Jiang [55]. Main features are summarized as following, 


1. General two-dimensional body-fitted coordinates. 

2. Time accurate transient capability. 

3. Conjugate Gradient Squared matrix solver. 

4. All speed flow capability. 

5. Pressure based method. 

6. k — e two equation turbulence model with wall functions. 

7. Non-staggered grid arrangement. 

8. Third order pressure damping term for face velocity calculation. 

9. TV D scheme with damping term. 

10. Cartesian velocity components solved. 


III.5.3 Two-Way Coupling Scheme 


The present study extended M — PI SO algorithm with more consistent updating in 
density variable and employed the concept of weak form transport equation [121]. 
The new procedure converged faster than our previous one [56] for compressible 
flows. A two-way coupling scheme is incorporated and presented in the following. 

The implicit finite difference form of the gas phase governing equations described 
in Chapter II can be written as: 
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i( /5 ” +i - n + a, ( P u, )> +i = sr/ 


(III. 44) 


- wri = n{vr') - A 'P” + 1 + s ? 1 + s:+j (m-«) 


■j pn- 1-1 pn 

±V.pH,)"”-(pH,r\ = G(tf,” +, ) + ^ + S^'+S^'+R}i'H^ (III. 46) 

^KA)” +i - (#»*)“] = /f(*" +i ) + sr 1 + (Hi-47) 
^[(/*) n+1 - (**)“] = £(t" +1 ) + S, n+1 + Sy 1 (III.48) 

^[(/>i/)” +1 - (pi/)”] = I(Y? +l ) + S3 1 - RT,' (HI-49) 

and 

- (pY„)'} = J(Y. 3 1 ) - 3 - (III. 50) 

where, the operators H, G, K, L, I and J stand for the finite-difference representations 


of the spatial convective and diffusive fluxes of momentum Ui, total enthalpy H t , 


turbulence kinetic energy k , its dissipation rate e, species mass fraction Yj and Y 02 
respectively, and the operator A, is the finite- difference equivalent of d/dx{. The 
droplet source terms and combustion source terms are listed in section (II. 4) and 
(III. 4), and the splitting of momentum source term in section (III. 3). 


The operator H can be constructed from C-0 TVD schemes with damping term 
for the convection terms [55] and the central difference scheme for diffusion terms. 



Generally, it takes the form 
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nb 


H(U t ) = A m U iim - A p Ui 


m = 1 


(III. 51) 

where tn is a grid node ld e„ tlfier ^ ^ ^ ^ ^ ^ 

nodes involved in the formulation of the finite-difference representation of the spatial 
^xes. Suffix , is for the centra, point. To ensure Better accuracy and stabdity of the 
overai, schetne, the centra, („., diagona,, Cement A, of the operator H is separated 

The rest of the Cements are retained on the right-hand side where they are treated 
P icitly m corrector stage equations. A new operator H' is defined as 

nb 

W) = ■ ff ( £/i) + A ” Vi = E fln.52, 

m=l 

Thus, a general implicit discretized transport equation can be written as 




nb 


( ^r + A ’)*r = E A ->K +i + + 5, 


The link coefficients A m and A 


771=1 


they are related 


p are instructed so that all of them 


(III. 53) 


are positive and 


as 


nb 


- 5T + (G e n+1 - C" +I + C" +1 - 


m=l 


C +I ) (III.54) 

’ C " ’ C” 41 C," +1 are face mass fluxes divided by the cell volume 

at east, west, north and south side respectively. A weak form of the transport 

equatton suggested By Yang et a,. [121] is obtamed By addtng the continuity equation 
multiplied by the dependent variable, 


~[ 


p n+1 - p n 


At + ( C e C Z +1 + CZ +1 - c; +1 )] • $”+! = -5 m 


(III. 55) 
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to the transport equation (111.53). If A, is replaced by 

nb 

Ap = ^2 Am + 


(III. 56) 


( i + A^r=i:^ +1+i£ sf +S4 


the resultant weak form of the transport equation becomes: 

nb t„<h \n ( 111 . 57 ) 

71+1 J ;fv71+l 

This weak form elim, nates term in the coefficient of *„ hence ft and *, 

are decoupled. The weak form is very helpful for constructing pressure correction 
equation due to this decoupling. The possible reason is that the term W* - *e 
coefficient of V?+' in the non-weak form represents a non-linear relation and may re 
suit numerical instability. Our previous numerical experiences with variable density 
calculations also show that the stability is enhanced for scalar transport equations 
by setting term to p" in the coefficient of *, for steady state calculations. 

Present discretisation is based on the explicit treatment of convective mass fluxes 

in transport equation, 

Sim + A 1( ^ r r«] = (UI 58) 

This explicit treatment avoids frequent calculation of coefficients, results a simple 
pressure correction equation and is numerically efficient. More investigations should 
be done for convection term treatment especially in time accurate transient calcu- 
lation of compressible flow and combustion process. The first term in the above 

equation is discretized as 

d (pg) _ P n $ n+1 - p n - x ^ n (III. 59) 

dt At 

to keep equation (111.58) satisfy Galilean transformation when a constant velocity 
or temperature is added. 
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The weak form equation (III.53) of the explicit 


convective flux treatment be- 


comes: 


1 rib 


m= 1 


At 


The two-phase momentum coupling source terms are listed again, 


(III. 60) 


<?n + i = J_ I 

<^1/ * ^ y-n+l 1 , At 

' L * *714-1 


p=l 


(III. 61) 


1 NP 4 Ar 
ffn+1 _ y- $*PdN p 

p ‘ Av2^—ZT~ f( r r) 


(C +1 ) 3 


At 

p=i 


At 


p) (u ’ )n * ~ FbiAi )~^ + v r]} (in. 62) 

and noting here, S p is independent of directions. 


Predictor Step 

Momentum equation is solved implicitly in this step, using pressure, density and 
two-way coupling source term quantities evaluated at the previous time step. 




( AT + W ^ ^ ' W + (III.63) 

Droplet injection, wall interaction, evaporation, breakup and collision can be ac 
tivated. Sub-time scale St„ is used for droplet evaporation calculation. Due to 
the nonlinear relationship of mass transfer number B m and droplet temperature T d , 
droplet temperature T d is solved iteratively using 


IT z II = A '"< r - T r)N . L » (pDrsmn(1 + fl , +1) 

Ste ° ~ (IIL64 > 


The droplet sire change will be calculated following the implicit solution of equa- 
tion (III. 64) for T v+1 , 
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, _ [r , ? _ (g£r 5 ^ fa (1 + B - ) + ' n(1 + B ^^ (Hi-65) 

v r P 1 -\'p> pd l 

At the same time, the source terms S„., and S w are evaluated using equations 
(III. 42) and (III.43). Particle velocity is updated, 


v i = 1 T At 

1 ~ T n 

These values of U- and u* are used to evaluate r*,S; and R p 


(III. 66) 


First Corrector Step 

We calculate temperature T‘ from relation 


1 


h; = h*(t) + -u:u: 


(III. 67) 


where total enthalpy H, is obtained from energy equation. Density is updated usmg 
the equation of state, 


pn 

nT _ _ 

9 RT * 


We approximate continuity equation as 


£l z£. + a t ( P nT ur) + &i(pun = s m ,i 

At 


(III. 68) 


(III. 69) 


where 


P = 


i?T* 


/ « nT 

P = P ~ P 


p » _ pn 

RT’ 


The momentum equation is approximated by 




(£ll + A P )U'“ = H\L !■) - A.P- + St + 

t 


(111. 70) 

(111. 71) 


- s:i/" + jc 


(111.72) 


59 


By subtracting equation (III.63) from equation (III. 72), we have 


on- 1 


( Ai + ^ + s;w - un = — Af(P* - n _ (5 ; - s;)U- + - ** (m . 73) 

Define the short notation 


«n— 1 


~ rj + a p + s;) 


-1 


(III. 74) 


Velocity U~ will be 


Ur = u; - D '^ p ' - n - Di[(s; - s;)u: - 


Substitution of equations (111.70,111.71,111.75) 
correction equation, 


( R ;< - k )i (in. is) 

into equation (111,69) yields pressure 


u: 


l AtRT- + Ai( #r) - d;a,)](p- _ = + a^u-)] 

+ V, + A. { ^; ((5; _ w _ (fl . _^ )J} (ni76) 

The particle velocity at this level is then 


v? + (Ur + u' + F bl r~)^ 

' l — * T * 


(III. 77) 


The values obtained at this level are used to calculate r~, 5;-. Al this stage , 
the mean velocity field satisfies the continuity constraint. 

Second Corrector Step 

To further sattsfy the momentum conservation, a second corrector step is used. We 
again calculate temperature T** from the relation 


H r = h"(t ) + l ~uru; 


(III. 78) 
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Density is updated using the equation of state, 


p mT = 


RT * 


(III. 79) 


We approximate continuity equation as 


{£ zZl + a t ( P - T ur) + Mpur) = 5m,/ 

At 


where 


p ** 

p " = 


*T _ 


p>** JD* 


P - P p ~ RT* 

The momentum equation is approximated by 


n n- 1 jyn 


ri — 1 P £»**Tf*** , 7 ">** 

(Ai + A P )t/r = ^(wr) - A * p " + Si + — “ s - • r ' 

By subtracting equation (III.72) from equation (III.83), 


(111. 80) 

(111. 81) 

(111. 82) 

(111. 83) 


„n-l 


{ ? + a, + s;-)(vr-un 

v A t 


= H\U** -U’)- A t (P** - P’) 

- (sr - s;)L r r + i?;; - ( IIL84 ) 


Define the short notation 


^ = ( ^r +Ap+5rrl 


(III. 85) 


Velocity I/*** will be 

y... = ur _ d”a,(p-- - n - cri(sr - s;)ur - (r; : - «;.)] ( 1IL86 > 

Substitution of equations (III.81, III.82, 111.86) into equation (111.80) yields pressure 


correction equation, 


1 


; — + A,(2E T )-A,(p- T i)rA,)](P"-p-) = -i^-^r L + A '('’' Tc, ''' )1 

l A tRT** RT- 

+ s m , - A ,[p- T D"H'[ur - v;)} 

+ A,{p- T i>ri(5r - s;wr - w - R ',.m 


(III. 87) 
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At this stage, the mean velocity field satisfies the continuity constraint. 

Following [47], it can be shown that the errors introduced by the operator- split- 
ting procedure is less than the truncation errors of the finite difference scheme used 
in the governing equations (III.45) and (111.15). Note that the effective relaxation 
tune r depends on the drag function which should contain the effects of turbulence. 
We therefore calculate (u t ) at this stage by a stochastic method with the k - c 
turbulence model. An one-predictor (im P licit)/one-correction (explicit) procedure 
for k and e equations suggested by Issa [47] has been used in this study. We then let 

V-" and P- be the value at *"+■ level and add the W) - to update the final time 
level particle velocities using the equation : 

^ + + + ]4L 

' ! ~ l + Ai (III. 88) 

' T* m 

The particle location is updated, 


= *1 + v? +1 At 


(III. 89) 


This brings all variables to the new time level. The time is then incremented and 
the new predictor-corrector procedure is repeated with the new velocities. This 
algorithm is used for simulation of transient phenomena. If only a statistically 
steady solution is desired, then the time steps for gas phase and particle phase can 
be made unequal j also the corrector stages of u, (III. 77 and III.8S) may be neglected. 

Scalar Equations 

After continuity equation 


— (p ^ 1 _ p«) + a i( P Ui) n+i = szy 


(III. 90) 
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is satisfied, we use the following weak form equation for scalars, such as energy, 
spec.es mass fraction and turbulent kinetic energy and its dissipation rate, 


r 1 *- 1 + 

At OXj t,x < 1 


Its discretized weak form is 


nb 




where the central point coefficient 


771=1 


nb 


Ap — ^ ^ Am “h Sm, l 


771 = 1 


The k and e discrete equations in the weak form are 

(El + C p )k n+1 = I<\k n+1 ) + + OtGT +1 - ( P e T +l + 5fc ’ i 

v A t ^ 


(III. 91) 


(III.92) 


(III. 93) 


(III. 94) 


and 


+ + D r )?*' = ty 


.71+1 


) + 


p n e n 

~KT 


+ (CiP(Ge/fc) n+1 - (Ciet'lk)"*' + S,j (11195) 


where 

nb 

r(fc" +, ) = ^c„c +1 

m= 1 
nb 

£'(e" +1 ) = E D ” e “ ‘ 

771=1 

ni> 

Cp = 'y ^ Cm H" 

771 = 1 
7lfr 

D p = + ■S’m.i 

m=l 

By invoking the eddy viscosity formulation: 

k 2 

Ht = C>— 


(111. 96) 

(111. 97) 

(111. 98) 

(111. 99) 


(III.100) 
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the source terms S k and £ are transformed into the forms 


and 


S k = {fi t G ) n+1 - (^k 2 ) n+1 

f*t 


s< = (c l c fi Gpk) n+i - (c 2 cy-e) n+i 


(III.101) 


(III. 102) 


Substitution of those changes, k and e equations become 

(£ + = *'(*«') + ££ + + S W 


(III. 103) 


and 


,k n+1 


{ At + Dr + C2C “ 1 ’ 2 '^ )e " +1 = + (CiC u Gpk) n+1 + S, j (III. 104) 

Predictor Step 

Decoupling the above equation as 

( A ~t +CpJr = K '( k ~) + + tfG + S kJ 

and 


(III. 105) 


( A t +Dp + C2C » p2 ^ € ' = L '( e "> + E XT + C ' c »Gp n+1 k' + S t j (III. 106) 

Equations (III. 103) and (III. 104) are to be solved in that sequence and implicitly in 
k * and e* respectively. Turbulent viscosity p t is then calculated from 

= C fi p n+1 tL- (III. 107) 

Corrector Step 

The corrector equations for k and e are explicit and can be written as 

( At + Cp + = *'(**) + + Pt ° + Sk ' 1 (HI. 108) 
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and 


(£ + D p + C 2 C»p 2 —)t' = 2/(0 + ^r- + C,C»Gp n+l k” + S e ,/ (III-109) 

At Pt 

Subtracting equations (III.108 and III.109) from equations (III.105 and III.106) re- 
spectively, we have 


p " , /-r , P O , P 


= <t; +c ’ + * 


*T'l(3j + c,+ 


P" 


k n )k* + (/q - /*?)G] (III.110) 


and 


= <X7 + D " + It + ^ + C2C *' ,! u? )e 


‘At r Pt 

+ C,C.Gp n +\k”-k')\ 


(III.lll) 


The new turbulent viscosity is obtained by 


Pt 


C.P 


n+1 


e 


** 


III. 5. 4 Summary of Solution Procedure 


(III.112) 


The numerical procedure implemented in the present MAST code is summarized as 


following, 


• Predict gas-phase velocities. 

• Calculate particle-phase properties. 

1. Inject particles. 

2. Find particle index with respect to grid points. 

3. Reflect particles if they cross the boundaries. 

4. Repeal particles which are out of calculation domain. 
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5. Calculate particle turbulent dispersion. 

6. Calculate droplet breakup and collision. 

7. Calculate droplet evaporation. 

8. Calculate two-way coupling source terms. 

• Correct momentum and temperature (two-steps). 

• Update particle velocities and trajectories. 

• Calculate species mass fractions. 

• Calculate total enthalpy. 

• Predict and correct turbulent kinetic energy and its dissipation rate. 

• Modify boundary conditions. 

• Go to the next time-step calculation. 

III. 6 Boundary Conditions 

For two-phase flows, both boundary conditions for gas-phase and droplet-phase must 
be specified to complete the solutions. 

HI. 6.1 Gas Phase Boundary Conditions 

Inlet boundary condition 

All dependent variables need to be specified at the inlet boundary. For incom- 
pressible flows, only pressure difference is required in calculations. Inlet pressure 
is obtained by linear extrapolation from inner point assuming same gradient. For 
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compressible flows, we need to specify pressure or stagnation pressure depending on 
subsonic or supersonic inlet. Detailed discussions can be found in Ref. [55]. 

Outlet boundary condition 

All of the variables are extrapolated at outlet boundary except velocity components 
in incompressible flow. Outlet velocity components are corrected based on the law of 
mass conservation, and inlet mass and droplet evaporation mass for incompressible 

flows. 

Symmetric boundary condition 

There is no flow across the symmetric line or axisymmetric axis and the gradients 
of solution variables are specified to be zero. 

Wall boundary condition 

No-slip condition is applied at impermeable walls. All gas velocity component and 
species gradients are specified to be zero. Temperature or heat flux can be imple- 
mented at walls for energy equation. For turbulent flow, almost all of the solution 
variables change sharply at near wall region. There are two alternatives to treat this 
stiff boundary condition. One is to use very fine grid in the laminar sublayer and 
buffer zone regions resulting in massive increase for computation time and computer 
memory, this method can be found in the low Reynolds number model of Jones and 
Launder [58] . The efficient method is the wall function treatment of Launder and 
Spalding [70], which is employed in the present study. 
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Open boundary condition 

Only pressure is specified at open boundary. The other variables are obtained from 
linear extrapolation. 

III.6.2 Particle Phase Boundary Conditions 

When particles are injected from a nozzle or atomizer, particle radius, temperature, 
density, velocity and mass flow rate need to be specified by experimental data, empir- 
ical correlations or testing data according to its operating conditions. These initial 
conditions are very important for spray dynamics, evaporation and combustion sim- 
ulations. A lot of experimental studies have been conducted for various atomizers 
to establish empirical correlations for spray flow rate, cone angle, size distribution, 
and penetration length [72], which are effected by liquid properties, atomizer and 
running condition etc.. 

Injecting particle size 

For dilute polydispersed sprays, particle radii are chosen from particle size distri- 
bution correlations, such as y - squard, Nukiyama-Tanasawa, or Rosin-Rammler 
distributions as described in Chapter II of Ref. [72]. If the injected computational 
particles and particle Sauter mean diameter, SMD or the related parameters are 
specified, we use the Monte Carlo method as described in section (III. 2) to choose 
each computational particle size. Calculation of each physical particle in a practical 
spray is impossible for current computer power. Based on the present statistical 
model, each computational particle represents a group of physical particles which 
have the same characteristic values, such as radius, temperature and velocity. Only 
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a limited computational particles are traced and calculated. Usually there are two 
kinds of method to represent physical particles if particle mass flow rate, size and 
computational particle number are specified. The first one specifies the physical par- 
ticle number jV p represented by each computational particle, and each computational 
particle has different mass 

m cp = | K{r p ) 3 pdN P (III. 113) 

Then, r p is generated using the drop number distribution function. In the second 
method, each computational particle has the same mass and represents different 
physical particle number. The mass can be obtained as 


FLOWP ■ At 
mcp ~ NPTS 


(III. 114) 


where FLOW P is the particle mass flow rate and NPTS is the computational parti- 
cle number injected per time step. Then physical particle number can be calculated 
with the following relation, 


1 


(III. 115) 


m 


cp 


and r p is generated from the drop volume distribution function. The second method 
is adopted in this study. 

For dense spray case, droplet breakup has to be considered. In the present 
study, both TAB model of O’Rourke and Amsden [89] and Reitz’s [95] breakup 
model are incorporated. In these models, the injecting particle sizes are the same 
as nozzle characteristical size. Droplet breakup occurrence depends on slip velocity, 
oscillation velocity and Weber number. 

In case of the droplet passage through the plane of symmetry, the mirroi image of 
the droplet with the same instantaneous properties, physical dimensions and velocity 
vector, is injected into the flow field. 



Droplet Wall Impingement 
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The wall impingement model adopts the jet treatment^] and the empirical corre- 
lation approach [84], The experimental data of Wachters and Westerling [118] can 
be numerically fitted in terms of the droplet Weber number before and after impact. 


^’ e ° ~ 0-678 We, exp(— 0.004415 We,) 


(III. 116) 


For We, < 80.0, the drops do not disintegrate during impact and bounce from the 
surface while for We, > 80.0, the disintegration produces a dispersion of the small 
drops on the surface. Thus, in case of We, > 80.0, the jet model is used; in case of 

< 80.0, the drops bounce from the surface and the normal velocity after impact 
can be calculated from the following equation, 


( 11 ) 


K = Vi(We 0 /W e,) 0 ' 5 

This wall impingement model is based on several assumptions such as extrapolation 
of the results with water drops at atmospheric conditions and at higher wall temper- 
ature, no breakup at impact, the neglect of the wall heat transfer, and the neglect 
of droplet interaction w.th a possible liquid wall film. Despite these limitations, 

the qualitative agreement for We, < 80.0 and the good quantitative agreement for 
We, > 80.0 have been reported. 

Injecting Particle velocity 


Ordinary atomisers produce a spray which is distributed in some cone angles with 0, 
and Let 0, and 0, represent inner cone angle and outer angle respectively, so we 
have 0, < 0,. If 0, = 0, particles are distributed throughout its volume and we call 
th,s kind of spray as sohd - cone spray. Otherwise, if particles are concentrated at 
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, h e outer edge of a conical spray pattern with * * 0 we call it hollo, - cone spray. 
If particle injection velocity V,„, and cone angles are specified from experimental 


data or empirical correlations and no more 


information can be used, we assume 


particles are uniformly distributed in the given 


cone 


angles. Axial velocity V x and 


radial velocity V. are chosen randomly from following expression, 


9 = 0 \ + ($2 $i)^/ 

V x = V in jCOs6 
V r = V in] sind 


(III. 117) 


where R , is a uniform distribution random function. For a swirling atomiser, particle 
velocity V z at tangential direction also need to be specified. 
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Results and Discussions 


In this chapter several benchmark problems involving particle dispersions in homo- 
geneous and inhomogeneous turbulent flows are studied to calibrate the stochastic 
method for particle-turbulence interactions, and to evaluate the present parcel PDF 
model. Non-evaporating transient and steady solid-cone sprays are used to test the 
models of dense spray effects including droplet breakup and collision. Turbulence 
modulation effect and the efficiency of parcel PDF model for a dense spray case are 
tested at a steady solid-cone spray. A transient hollow-cone spray shows compli- 
cated two-phase interactions. Transient evaporating and burning sprays incorporat- 
ing droplet bieakup model, evaporation model and eddy- breakup combustion model 
demonstrates the sophisticated structures of such polydispersed spray combustion. 
The predictions show reasonably agreements with available experimental results. 

IV. 1 Particle Turbulent Dispersion 

IV.1.1 Nearly-Homogeneous Turbulent Dispersion 

The particle dispersion experimental setup of Snyder and Lumley [113] in a grid- 
generated turbulent flow was used for evaluating the present particle dispersion 
model and parcel PDF model. The mean flow was uniform in the test region and 
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turbulence was also found to be nearly isotropic. Individual particles were isokinet- 
ically injected by air flow through the sampling tube and carried to the center of 
the tunnel. The particle concentration was so low that mean gas flow properties 
were not effected. The turbulence energy decay curves were presented as a function 
of axial position. In the present study, turbulent kinetic energy was obtained by 
averaging axial and transverse fluctuating velocity and was fitted to the measured 
correlation as: 


k(x) = 


1.5C/ 2 


(IV.l) 


54.88 • - 14) 

where U = 6.55m/s is the mean axial velocity, x is axial distance from the mesh, and 
M = 0.0254m is the mesh size. The dissipation rate t is calculated by differentation 
of the above relationship, 


dk 

e(x) = -U- = 


l.oU 3 


(IV .2) 


dx 54.88M • (fj — 14) 2 
Particle pictures were taken at ten stations in the experiment, spaced logarithmically 
from x/M = 68.4 to 168. 

Generally, particle motion in turbulent flow is governed by the coupled effects 
of the turbulent flow field, particle inertia and crossing trajectory. All the three 
effects could happen for the dispersion of heavy particles, while the dispersion of 
light particle is dominated by the effect of turbulent flow field where light particle is 
fully correlated with turbulent fluctuating. For the motion of heavy particles, due to 
the influence of the Earth’s gravitational field, a particle’s free fall velocity increases 
with its inertia, which results the crossing trajectory effects. So, both effects of 
inertia and crossing trajectory are coupled together and hard to separate. In this 
experiment, particle densities and sizes are chosen to examine the three effects. For 
the light particle (hollow glass) with the effect of turbulent flow, the eddy lifetime 
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Table IV. 1: Particle Parameters 



Parameter Hollow glass Solid glass Corn Copper 


Diameter (p) 


46.5 


87.0 87.0 46.5 



Density (g/ cc ) 


0.26 


2.5 1.0 8.9 



controls the interaction times . For the medium particle (solid glass), the part.cle 
inertia or the transit tune controls the interaction times. But lor the heavy particles 
(corn pollen and copper), the crossing trajectory effect dominates flu. dtspersmn. 

this experiment, the concentration of petiole is so dilute that only one-way 
coupling is of main concern. In the calcnlation, the experimentally determined tad 

turbulence kinetic energy and its dissipation rate are used. The eddy hfe time 

i • o +• /tt and the eddy length scale l e = 1-65CV t 

t _ 3 was derived in Section (ll.d.4j ana rne j 

„ , • r al j ata The particle calculations were started at 

was tuned to fit the experimental data, in P 

the experiments particle injection point of ,/* - 20 . The part.de velocities 
were assumed e q ua, to the mean tad veiocity of 6,5 m/s and their fluctuating 
velocities were set to aero. When the part.cles were injected, each particle mteract.on 
time with eddy was randomly assigned within the eddy life time. Otherw.se the 
predicted particle dispers.on curve would be oscllatory with the period of eddy hfe 
time especially for light particle (hollow glass). 
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For the delta function SSF computations, 5,000 computational particles were 
sampled to calculate the resulting mean squared dispersion with respect to time. 
This computational particle number was chosen to get a smooth dispersion curve 
with minimum required number. 

For the PDF computations, a single parcel in a deterministic trajectory along 
the centerline was sampled to evaluate the mean squared dispersion representing the 
variance of the parcel PDF by using the related parameters for each eddy interaction. 
The calculation of this single particle was started with x/M — 68.4, where particle 
pictures were taken. The mean particle injection velocity was same as the fluid 
velocity, but its root mean squared ( rms ) velocity was set to 60% of the gas rms 
velocity to partially take into account the turbulence effects from x/M = 20 where it 
was injected. 

Figure IV. 1 shows the comparison of the predicted and measured particle dis- 
persion with respect to time. The particle dispersion is defined as the mean squared 
displacement y related to the centerline, 

. NP 

y! = i (iv.3) 

n=l 

The parcel PDF results show good agreement with the SSF results for light, 
medium, and heavy particles. Both models also show favorable agreement with the 
experimental data. 

These numerical results indicate that the parcel PDF model can accurately and 
efficiently predict the particle dispersion this nearly-homogeneous turbulent flow. 

IV. 1.2 Inhomogeneous Turbulent Dispersion 
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Figure IV. 2: Grid system and boundary conditions for round turbulent jet. 
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Figure IV. 4: Normalized particle concentration distribution ot P article la £ eD ™ 
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Figure IV. 6: Normalized particle concentration distribution ol particle laden round 
jet for PDF (200 parcels) with various correction factors (x/d— 20). 





Figure IV.7: Normalized particle concentration distribution of particle laden round 
je tor PDF (200 parcels) with various correction factors (x/d=30). 



Figure IV. 8: Normalized particle concentration distribution of particle laden round 
jet for PDF (200 parcels) with various correction factors (x/d=40). 






83 


The next example problem is the particle laden round jet of Yuu et al. [123] in which 
the turbulence is inherently inhomogeneous. The air jet is directed upward against 
the gravitational force. The nozzle diameter is 8 mm with a mean exit velocity 


of 20 m/s. Particles are injected with air velocity from the nozzle. The particle 
concentrations at the nozzle exit were measured at the centerline in the potential core 
and were between 0.8 and 4 ,/m’. The loading effects can be neglected, hence only 
one-way coupling needs to calculate. The turbulent gas-phase transport properties 
are provided by using the * - e model. A 41x31 grid system with uniform grids in 
axial direction and clustered grids about the centerline is used in the calculation. 

g e IV.. shows the vicinity grid distribution in the centerline region and boundary 
conditions. The injected turbulent kinetic energy 1.2, which corresponding to 0.3% 
U 2 , and its dissipation rate 20,000 are used to keep a very low eddy viscosity. Uniform 
profiles of velocity, turbulent kinetic energy and its dissipation rate are assumed for 
inlet boundary conditions. Numerical experimental results show that this choice 
gives a good prediction of the potential core length of the air jet when compared 
with the experimental data. For particle tracking calculation, a constant time-step 


lO/is is used. 


The particle concentration profiles have been extensively examined by Shuen 
[107] using the SSF model. Also, Litchford and Jeng [75] have tested their DDWT 
and SDWT models with similar flow properties and different particle sizes, and the 
SDWT model is suggested to take into account the flow property variation in radial 
direction. The present study is to incorporated the parcel PDF (or SDWT) model 
and compare with the SSF model in terms of accuracy and efficiency in the statistical 
model, and no intention is made to compare with the experimental data. In the 
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following calculation, the particle size is 20 fim and density is 2.0 5/cm 3 . Particle 
number concentration profiles are obtained with the particle number located in the 
cell divided by its volume and normalized by the total particle number m the cross 

section. 

Figures IV. 3, IV.4, IV. 5, IV. 6, IV. 7 and IV. 8 show the particle concentration 
profiles of the delta function SSF model and the SDWT model for 50 and 200 com- 
putational parcels, at various levels of the correction factor, and at several axial 
locations. 10,000 particles are sampled for the delta function SSF computations. 
Using the 10,000 particles in the SSF model, there is still evidence of slight under- 
sampling. However the distribution is relatively smooth and is taken here as a good 
approximation to the theoretical profile. The results of the SDWT model with 50 
parcels shown in Figures IV.3,IV.4 and IV.5 are very sensitive to the level of the 
correction factor, especially for upstream regions due to undersampling. By increas- 
ing the correction factor, K in equation (11.89), the uncertainty level m the mean 
increases the dispersion and smoothes the profile considerably. In Figures IV.6, IV.7 
and IV. 8 the zero correction factor case(K=0) corresponds to the delta function SSF 
case using 200 computational parcels. The computed profile of the SSF model using 
200 parcel samples is very irregular and shows oscillatory distribution. 

The 200 parcel case of parcel PDF model shown in Figures IV.9, IV.10 and IV.ll 
is less sensitive to the correction factor since there is less uncertainty in the mean 
because of increased sampling. In Figures IV.9, IV.10 and IV.ll, the PDF results 
with 200 parcels and K=4 show favorable agreement with the delta function SSF 
with 10,000 computational particles. In terms of the CRAY X/MP-24 CPU time, 
the parcel PDF (SDWT) solutions with 200 parcels requires about 36 seconds while 
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the SSF solutions with 10,000 parcels needs about 1,375 seconds. These numerical 
results clearly indicate that the parcel PDF model has the capability of accurately 

representing dispersion in inhomogeneous turbulent flows with improved efficiency 
over the delta function SSF model. 

In the calculation of Litchford and Jeng [75], with only 50 parcels the PDF model 
produces smooth and accurate distributions comparing 200 parcels required in the 
present calculation. This discrepancy should be attributed to the different methods 
for calculating the particle concentration and different grid size used for accounting 
the particle number in the cell. In Litchford and Jeng’s calculation, they used 
PSIC and particle concentration was obtained by calculating the particle transition 
time crossing the cell boundaries, which usually gives more accurate results but 
needs more computation time and is not convenient for arbitrary grid system. In 

the present method, only the particles which located in the cell are taken into the 
concentration calculation. 

IV.2 Non-Evaporating Solid-Cone Spray 

IV.2.1 Measurements of Hiroyasu and Kadota 

The solid-cone spray measurements of Hiroyasu and Kadota [43] were used to 
validate the present numerical dense spray model which includes collision, coales- 
cence, and breakup models described in Chapter II. Liquid fuel is injected through a 
single hole nozzle into constant pressure, room-temperature nitrogen. Spray tip pen- 
etration and drop sizes were measured from photographs of the backlighted spray. 
The test conditions are given in Table IV.2, and the SMD is measured for the spray 
cross-section 65 mm downstream of the nozzle. The nozzle diameter was 0.3 mm 
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Table IV. 2: Test Conditions for the 
Nozzle diameter: 300 yum 

Liquid: diesel fuel 

Viscosity: 3.8 x 10 _4 jV • s/m 2 
Ambient gas: nitrogen 


Measurement of Hiroyasu and Kadota 
Injection pressure: 9.9 MPa 

Density: 840 kg/m 3 

Surface tension: 0.0232N/m 
Temperature: 25 °C 


Case 


Pgas 

(MPa) 


Pgas 

(kg/m 3 ) 


Vinj 

(m/s) 


Minj 

( k g/s) 


SMD 

(/zm) 



1 

2 

3 


1.1 

12.36 

115.80 

0.00688 

42.4 

3.0 

33.70 

102.54 

0.00609 

49.0 

5.0 

56.17 

86.41 

0.00513 

58.8 



“ d Pr ' Sent tetrad ecane for the ,i q uid fee, (the experiments 

used a diesel fuel with physical properties close to tetradecane). 

A computational domain of 20 mm in radius and 120 mm in length was discretized 
al and 45 axial grids. The mesh spacing was nonuniform with refinement 
on the centerline and close to the injector. The smallest cell is 0.5 mm radially and 
1.5 mm axially. The grid system and boundary conditions are shown on Figure 
IV.12. The time-step sizes used for three test cases are 10. 18 and 18 ps respectively 
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and boundary conditions for solid-cone spray. 
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and 250 time steps are used for all of the three cases. The number of computational 
parcels at the end of calculations was between 1000 and 3000, which were varied 
with the back pressure. The present numerical results did not change appreciably 
when this parcel number was varied. The initial turbulent quantities were assumed 
as the small values (k = 1 x I0~ 3 m 2 /s\t = 4 x 10- 4 m 2 /s 3 ). The numerical results 
were insensitive to these initial values. The typical CPU time required on CRAY 
X-MP/24 using both breakup models is listed on Table IV. 3. 


Table IV. 3: CPU Time Requirement for the Hiroyasu’s Case 



TAB Model 

Reitz’ 

s Model 

Case 

Parcel 

CPU Time (s) 

Parcel 

CPU Time (s) 

i 

1130 

114.6 

1854 

207.0 

2 

1181 

120.4 

1038 

150.4 

3 

1248 

123.4 

1676 

229.0 


Figures IV.13, IV. 14 and IV.15 show the spray parcel distributions for three cases 
with the TAB breakup model [89] and Figures IV.16, IV. 17 and IV. IS with Reitz’s 
model [95]. These plots indicate that the spray tip penetration and the core length 
decrease with the increase of the gas density. The predicted core lengths (20mm for 
P=l.l MPa, 12mm for P=3.0 MPa, 10mm for P=5.0 MPa) based on Reitz’s model 








94 



Figure IV.15: Spray parcel distribution in a solid-cone spray (Case 3, Time-4.5ms, TAB model). 
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have a reasonable agreement with the jet length correlation suggested by Arai et 
al.{6], but the TAB model predicts much shorter intact core lengths. These results 
indicate that the Reitz’s wave instability model is capable of predicting the intact 
core length with re£ison3-ble ctccur&cv. 

Figures IV. 19 and IV. 20 show the predicted and measured spray tip penetration 
for two breakup models. The calculated penetrations for two breakup models show 
reasonably good agreement with the measurement. Compared to the TAB model, 
the Reitz’s model slightly overpredicts the penetration length for three cases. How- 
ever, the discrepancies in the penetration length could be partially attributed to the 
imprecise definition of the spray tip. In the present computations, the spray tip was 
defined to be the location of the leading spray drop parcel. It is necessary to note 
that a far-field spray penetration is not a sensitive indicator of model performance. 
Previous studies [31, 21] indicated that a far-field spray penetration is mostly influ- 
enced by the turbulence diffusivity. However, a near-field spray penetration could 
be more sensitive to the physical submodels such as breakup and collision. 

Figures IV.21 and IV.22 show the variation of SMD for two breakup models. 
The three data at 65 mm correspond to the measurements. The computed drop 
sizes are time-averaged over the spray cross-section at each axial location. At the 
nozzle exit, the drop diameter is equal to the nozzle diameter, 0.3 mm. The overall 
trend of the SMD distribution is similarly predicted by the two models. Generally, 
these curves can be broken into two sections. Close to the injector, the drop size 
decreases rapidly due to drop breakup. Further downstream, the drop size increases 
gradually due to drop coalescence. In the lower pressure case(l.l MPa), the drop size 
remains relatively uniform after initial breakup region and then increases slightly in 
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IV .20: Spray tip penetration 




103 


the far-downstream region. For the high pressure cases(3.0 and 5.0 MPa), the drop 

S1Z e largCly m downstream regton, because higher gas densities promote 

collisions and coalescence. This trend is also observed in the measurements. The 
predicted drop sizes at 65 mm qualitatively agreed with the experimental data for all 
three cases. The discrepancy could be associated with the fact that the experimental 
sprays were pulsed while the computations assumed a constant pressure injection 
for the entire computational time period. In comparison with the Reitz’s model, 
the TAB model provides the relatively rapid breakup rate near the injector and 
the relatively low SMD distribution for three cases. Especially close to the injector 
for the low gas pressure casefl.l MPa), the TAB model predicts the much faster 
breakup than the Reitz’s model. Due to this faster breakup rate near the injector, 
the TAB model underpredicts the intact core length for three cases. The larger SMD 

distribution predicted by the Reitz’s model could be tied with to the overprediction 
of the spray penetration length. 

IV.2.2 Measurements of Wu et al. 


To compare with the local flow properties in terms of the velocities and the 
turbulent intensities for gas and droplet, the measurements of Wu et al. [120] were 
selected. In this experiment, axial and radial components of the droplet velocity 
were measured by laser Doppler velocimetery (LDV) within liquid n-hexane sprays 
injecting into high-pressure nitrogen from single-hole cylindrical nozzles at room 
temperature. It was found that beyond 300 nozzle diameters from the nozzle that 

the fully developed incompressible jet structure and droplet-gas equilibrium were 
being approached. 

Because this is a steady state case, no initial conditions are required. The grid 
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Table IV.4: Test Conditions for the Measurement of Wu et al. 
Nozzle diameter: 127 Injection pressure: 9.9 MPa 

Liquid: n-hexane Density: 665 kg/m 

Viscosity: 3.2 x 10~ 4 N ■ s/m 2 Surface tension: 0.0184 N/m 

Ambient gas: nitrogen Temperature: 25 °C 


Case 

Pinj 

(MPa) 

Pgas 

(MPa) 

Pgas 

(kg/m 3 ) 

Vinj 

(m/s) 

x/ d 0 

A 

12.5 

1.48 

17.02 

127.0 

400 

B 

15.2 

4.24 

48.68 

127.0 

500 

C 

30.4 

4.24 

48.68 

194.0 

500 


system and boundary conditions are the same as the previous case, shown on Figure 
IV. 12. The time-step sizes used for three test cases are 20, 20 and 10 \is respectively. 
At the beginning of the calculation, larger time-step sizes are used to establish the gas 
flow field, and then the above small time-steps are used for two-phase calculations. 
300 time steps are used to obtain statistically steady results. The droplet properties 

and three injection conditions are listed on Table IV.4. 

In this experimental setup, the gas turbulent round jet was developed by the 
injection of liquid fuel and its atomization. Two-phase flow was fully coupled by 
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gas-droplet momentum exchange especially in the immediate vicinity of the nozzle 

exit. For the present dense spray, it is of interest to establish the importance of 

droplet modulation to gas turbulence and compare with the available experimental 

measurements. Before getting into the detailed comparison with the experimental 

effect of drops on the turbulence is studied with several modulation models. 

In the following figures, Model I corresponds to the model suggested bv Mostafa 

and Mongia [83]; Model 2 by Amsden et al. [4] and Model 3 by Chen and Wood 
[78], 


Figure IV 7 . 23 shows 


the predicted gas-phase centerline velocity with and 


without 


modulation terms in standard k - e equation. Close to the nozzle, all models predict 
the peak velocity at the nearly same axial location. However, there exists differences 


Of the centerline velocity distribution near the 


nozzle exit. At downstream of the 


peak point, the standard * - . mode, without modulation effect predicts the fastest 
decay of the centerline velocity followed by Model 2, Model 1 and Model 3. In the 
far downstream, all models predict similar slopes due to the relatively small effects 
Of the turbulence modulation. Therefore, the turbulence modulation is important in 
the jet developing region. Model 1 and Model 3 show the favorable agreement with 
the far-field measured velocities(Case C). Figure IV.24 shows the centerline turbulent 
kinetic energy with four models. Around the peak point, Model 1 predicts the lowest 


turbulence level followed by Model 3, Model 2, and no-modulation case. In the far 
field, Model 1 and Model 3 have the almost same turbulence level. 

Figure IV.25 shows the comparison of predicted and measured centerline velocity 
for the three test conditions. Predictions based on Model 1 have a good agreement 
With measurements for three test cases. Figure IV.26 shows the radial profiles of gas 


Centerline Velocity (m/s) 
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Figure 


IV.23: Gas-phase axial centerline velocity (Case C) 


Centerline TKE (mV's 2 ) 


7 



Figure IV. 24: Gas-phase axial 


centerline turbulent kinetic 



Centeriine Velocity (m/s) 
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Figure IV.25: Gas-phase axial centerline velocity for three test conditions 
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Figure IV.26: Radial profiles of gas/drop 
Case A) 


axial mean and RMS velocity (x=50.8mm, 
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and drop axial mean velocity and RMS velocity at axial location of 50.8 mm (400 
nozzle diameters, Case A) from the injector. The predicted mean drop velocities 
weighted by the drop number are obtained by averaging the instantaneous values 
over 6ms period. The predicted mean and RMS velocities are somewhat lower than 
measured database A). As would be expected, the predictions for gas and drop 
velocities are in almost equilibrium at this downstream location. 

The predicted and measured results for Case B are shown in Figure IV .27. The 
corresponding axial location in Case B is 63.5 mm (500 nozzle diameters) from 
the injector. The computed mean and RMS velocities are favorably agreed with 
the experimental data. The slight oscillations in the predicted drop velocity profile 
are due to undersampling of computational parcels. Figure IV.2S summarizes the 
results for Case C at axial location of 63.5 mm from the injector. In this case, the 
overall agreement is good in the mean and RMS velocities. However, around radial 
location of 4.5 mm, the predicted mean drop velocities are somewhat lower than the 
measurements. The previous numerical study of Reitz and Diwakar[97] showed the 


similar trend with the present study. 

Figures IV .29, IV.30 and 1V.31 show the radial profiles of the mean gas and drop 
velocities, and the instantaneous drop velocities using the SSF model and the parcel 
PDF model for three cases. The number of computational parcels for the SSF model 
and the parcel PDF model is about 3600 in the whole flow field. Compared to the 
parcel PDF results, the computed profiles of the SSF model for three cases are very 
irregular and oscillatory. The parcel PDF model provides the realistic and similar 


distributions with the mean drop profiles. Due to slight undersampling, certain level 
of irregularities exist in the distributions of mean drop velocities and instantaneous 
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Figure IV.29: Radial profiles of 

comparisons with PDF (Case A) 
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Figure IV.30: Radial profiles of gas/drop mean and instantaneous axial velocity 

comparisons with PDF (Case B) 
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Figure IV.31: lUdial profiles of gas/d 

comparisons with PDF (Case C) 
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droplet velocities using the parcel PDF model. Compared to the SSF calculations 
with the same number of computational parcels, the CPU time with the parcel PDF 
model is increased about 10%. However, to reach the same level of the realistic 
distributions of the parcel PDF model, the SSF model needs to increase the num- 
ber of the computational parcels at least several times which results in the drastic 
increase in CPU time and the memory storage requirement. Furthermore, the irreg- 
ular drop/gas distributions due to the undersampling of the computational parcels m 
the SSF model can create numerically oriented high-frequency noises which possibly 
contaminate the high-frequency combustion instability solutions. The present parcel 
PDF method has potential advantage to eliminate numerically oriented noises asso- 
ciated with the droplet injection as well as to reduce the number of computational 
parcels with accurate representation of spray dynamics. The results indicated that 
the present parcel PDF model has the capability of accurately representing drop 
dispersion in dense sprays with manageable number of computational parcels. 

IV. 3 Non-Evaporating Hollow-Cone Spray 

The hollow-cone spray tip penetration data of Shearer and Groff [105] have been 
used for the model validation. In the experiment, the liquid fuel (indolene-clear 
gasoline) is injected into quiescent room-temperature nitrogen at P = 550 kPa. 
A computational domain of 20 mm in radius and 40 mm in axis was uniformly 
discretized by 41 and 81 grids respectively as shown on Figure IV.32 with boundary 
conditions. The experimental spray cone angle is 60 degrees, and the flow rate 0.0165 
ml/injection with four pulses, each of duration about 0.58 ms. The droplet injection 
velocity could be approximated as a constant through the injection [96] and is set 
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Table IV. 5: Test Conditions for the Measurement of Shearer and Groff 


Pgas 

(kPa) 

Pgas 

(kg/m 3 ) 

Vinj 

(m/s) 

VOLinj 

(ml/inj) 

Cone Angle 
(deg) 

550 

6.36 

60 

0.0165 

60 


to the experimental spray tip velocity (60m/s) measured from the movie pictures 
in the early stage of the injection. Turbulence modulation model and PDF model 
are not activated in this calculation. The TAB breakup model and collision model 
axe used to account for the effect of droplet breakup and collision/coalescence. 10 
droplets axe injected each time step with the same size as the nozzle drameter 60 
pm. The numerical timestep size is 10ps and about 1200 spray parcels are present 
at time 1.32 ms (132 time steps) in the computation with 141 second CPU time on 

CRAY X-MP/24. 

Figures IV.33,IV.34,IV.35,1V.36, IV.37 and IV.38 show the spray parcel distribu- 
tions and the velocity vectors with time advancing. The numerical results indicate 
that turbulence has a relatively small effect on penetration in a hollow-cone spray 
because the radial spreading due to inertia is the dominant factor. The gas velocity 
vectors indicate the presence of a vortex near the head of the spray, which curls the 
spray tip toward the outside of spray. A substantia, region of strong inward ffow 
in the center of the cone near the injector was also observed. These flow patterns 











Figure IV. 38: Velocity vectors in hollow-cone spray (Time= 1.32ms) 
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and spray shapes clearly show a very complex structure of a hollow-cone spray and 
two-phase interactions. These results are compared favorable agreement with the 
experimental observations [105]. The predicted and measured spray tip penetration 
versus time are shown in Figure IV. 39, and the present predictions show reasonable 
agreement with the experimental tip penetration. 

IV. 4 Evaporating and Burning Solid- Cone Spray 


Table IV. 6: Test Conditions for the Measurement of Yokota et al. 


Case 

Pinj 

Pgas 

Tamb 

Minj 

Atmosphere 


(MPa) 

(MPa) 

(K) 

(kg/s) 


Evaporating 






Spray 

Burning 

30 

3.0 

900 

0.00326 

N 2 

Spray 

30 

3.0 

900 

0.00326 

Air 


The evaporating and burning solid-cone spray measurements of Yokota et al. 
[122] have been used to validate the present numerical dense spray model with 
droplet evaporation. Liquid fuel (tridecane C 13 // 2 S ) is injected through a single hole 
nozzle into high-pressure, high-temperature nitrogen or air. The test conditions for 
evaporating and burning sprays are given in Table IV. 6. The nozzle diameter was 
0.16 mm. A computational domain of 20 mm in radius and 100 mm in length was 
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Figure IV. 39: Spray tip penetration versus time in a hollow-cone spray 
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Figure IV. 40: Grid system and boundary conditions of evaporating and burning cases 



Table IV. 7: Properties of tridecane 


Pd 

840 kg/m 3 

H corn 

4.45 xlO 7 J/kg 

MW 

184 kg/kmole 

T cr 

677 K 

T d 

300 K 

Vinj 

185 m/s 

a 350 

0.02126 N/m 




Table IV. 8: Constants for gas property 


AIRLA1 2.52 x 10~ 3 lV/(m • I<) AIRLA2 200 
AIRMU1 1.46 x 10" 6 ^/(m-s) AIRMU2 110 


discretized by a 21 radial and 44 axial grids. The mesh spacing was nonuniform with 
refinement on the centerline and close to the injector. The grid system and boundary 
condition are shown on Figure IV. 40. The number of computational parcels at 
the end of the calculations was about 400 and 900. More parcels have been tried 
with no significant effect on the calculation results. Due to the numerical reasons, 
the initial turbulent quantities were assumed to be small values (k = 0.2 m 2 /s 2 
and e = 1.0m 2 /s 3 ). The upstream boundary is treated as a solid wall, and other 
boundaries are treated as open boundaries. Collision model, turbulence modulation 
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model and PDF model are not activated in these calculates. A constant time-step 

"" 10 "* ““ 6 ^ "" ^ P- t^e-step for both evaporating 

and burning cases. After 4 ms of droplet injection with 400 time steps, the droplet 

numbers in the flow field were about 930 and 410, and the CPU time on CRAY X- 

MP/24 were 212 and 188 seconds for evaporating and burning sprays respectively. 

P perties of liquid fuel (tndecane C 13 // 26 ) are listed on Table IV. 7. Where 

case is the droplet surface tension at 350 K. The physical properties of surface tension 

fuel vapor diffus.vity multiplying by density P D, gas-phase thermal conductivity 

K, and gas-phase viscosity „ are calculated using the following empirical relation as 
functions of temperature, 


a = 


Ter ~ T d 
T cr - 350 a350 


(IV. 4) 


pD = AIRDIF ■ T expdif 

9 


(IV. 5) 


K = 


A1RLA\ ■ T'PKT, + A1RLA2) 


(IV.6) 


P = AIRMU1 - r; V(T s + AIRMU2) (IV 7) 

where T t and are droplet and gas-phase temperature respectively. The constants 

droplet and fuel vapor pressure as functions of temperature are taken in table from 
JANAF data bank [116], 

IV.4.1 Evaporating Solid-Cone Spray 



Figure IV.41: Spray parcel distribution in an evaporating spray (Time-0.2ms) 





















Figure IV.49: Contour of fuel mass fraction in an evaporating spray (Time 4.0ms) 
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Using the numerical results of the Reitz’s breakup model [95], the spray parcel 
distribution, the contours of the fuel mass fraction, and the temperature contours 
for the evaporating spray are plotted in Figures IV.41, IV.42, IV.43, IV. 44, IV.45, 
IV. 46, IV.47, IV.48 and IV .48 at different times of spray development. These results 
show that the spray penetrations increase with respect to time at early period of 
injection, however the penetration become nearly constant after t = 0.2ms due to 
evaporation. Even though the liquid drop does not penetrate more, the evaporated 
fuel vapor continuously penetrate with respect to time. At later period of injection, 
the low-temperature zone near the injector is created due to the cooling effects 
of evaporation. In this evaporating case, we did not turn on the collision routine 
because the collision process causes a significant increase of the drop size and the 
corresponding penetration length is unrealistically long. The discrepancy could be 
attributed to the inconsistency between the collision model and the vaporization 
model as well as the neglect of the supercritical vaporization effects. A similar spray 
characterization was obtained using the TAB breakup model. 

Comparisons of the computed and experimental spray penetration versus time 
are shown in Figure IV .50. Due to the neglect of the collision process, two breakup 
models underpredict the penetration length. At the initial stage of the injection, the 
TAB model noticeably underpredicts the penetration length. The underpredicted 
penetration length with the TAB model results from the rapid breakup rate near 
the injector. 

IV.4.2 Burning Solid-Cone Spray 

When the experimental atmosphere changed from nitrogen to air, ignition and com- 
bustion occurred at this high temperature environment. A single-step fast chemical 
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Figure IV. 49: Spray tip penetration versus time in an evaporating spray. 
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Figure IV.51: Whole net heat release rate of fuel vapor versus time. 
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reaction and eddy breakup model mentioned in Section (II. 2. 4) is employed here. 
The reaction for tridecane can be expressed as: 

C 13 H 28 T 2 OO 2 = I 3 CO 2 + IAH 2 O (IV. 8) 

The chemically controlled reaction rate for the fuel, given by the usual Arrhenius 
formula [119], can be written as following after units changed from cm- mole to m-kg 
system, 

- « w > ■ < IV - 9 > 

where the units of R c h. e is kg/(m 3 ■ s), p is kg/m 3 , E/R is Kelvin and molecular 
weight is kg/kmol. The constants used above are listed on Table IV. 9. Due to 
the conservation of atoms, only two mass fractions need to be solved. Fuel and 
oxygen mass fractions are chosen as dependent variables and the remaining ones are 
determined from the stoichiometric relations in Section (II.2.4). The stoichiometric 
constant s is calculated as: 

W 0 2 

5 = 20—^ = 3.48 (IV.10) 

Wf ' 


Table IV. 9: Constants in Arrhenius law 


Wf 184 kg/kmol A 3.2 xlO 11 a 0.25 

W 02 32 kg/kmol E/R 15,100 I\ b 1.50 


When liquid fuel jet is injected from nozzle, a spray is produced and droplets 
evaporate resulting fuel vapor accumulated. Ignition is accompanied following rich 
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vapor environment and high-temperature. The net heat release rate of fuel vapor in 
whole flow field is calculated and shown in Figure IV. 51. Before time=1.0 ms, the 
net heat release rate is very low and that means no combustion occurs. It increases 
dramatically after time=1.0 ms and reaches its maximum value which is much lager 
than the fuel supply rate. This should be attributed to the accumulated fuel vapor 
reacting with oxygen rapidly after the occurrence of ignition. The predicted ignition 
time 1.0 ms is in good agreement with the experimental results observed from the 
Schlieren photographs. The net heat release rate decreases from t=1.3 ms to 2.0 ms, 
that means accumulated vapor consumes a lot due to the rapid spreading of com- 
bustion wave. The increase of the net heat release rate after reaching its minimum 
value can be explained to the enhanced droplet evaporation rate by higher ambient 
temperature after the combustion. 

Figure IV. 52 compares the penetration lengths of evaporating and burning sprays. 
Before time=1.0 ms of ignition occurring, both penetration lengths are almost the 
same because the gas properties of nitrogen and air are quite similar. The burn- 
ing spray penetration length decreases after of ignition due to the higher ambient 
temperature promotes droplet evaporation rate. Both lengths keep almost uniform 
values after certain times, that means quasi-steady states are reached for liquid 
sprays even though the fuel vapors and temperatures still penetrate with respect to 
time. 

Figures IV.53, IV.54, IV.55, IV.56, IV.57, IV.58, IV.59, IV.60, IV.61, IV.62, 
IV. 63 and IV. 64 show the spray parcel distribution, the contours of fuel mass frac- 
tion, temperature and oxygen mass fraction at different times of injection for burning 
sprays. These figures clearly illustrate the relative locations of the cold core region, 
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reaction zone and spray boundary of a coaxial spray diffusion flame, and their devel- 
opment with time progressing. The liquid jet leaving the injector fully breakups very 
quickly and the spray is highly nonuniform with small droplet evaporating up after 
a short distance from the injector due to high temperature environment. Near the 
injector exit, the two-phase slip velocity is the greatest and the momentum of the 
spray is transferred to the gas over an extended axial distance. The small droplets 
around the periphery of the spray exchange the momentum with the gas and cause 
the spray jet to entrain surrounding gas, which also providing oxidizer for fuel vapor 
combustion. At the same time, the small droplets evaporate rapidly to provide fuel 
vapor which is consumed near the outer portion of the turbulent diffusion flame. 
The computed configuration of a spray flame has the overall agreement with the 
experimental observations. The comparisons of the predicted flame front movement 
with time for this transient spray combustion flow also show reasonable agreement 
with the experimental Schlieren photographs. 

In the experimental study, a considerable level of soot was observed near the 
spray tip where the equivalence ratio is low and the temperature is high due to the 
progressed turbulent mixing. Therefore, the soot model should be incorporated to 
improve the prediction capability of the present burning dense spray model. Future 
studies may include the detailed comparison with the local properties available in 
the experiment. 



V 


Conclusions and 
Recommendations 

V.l Summary 

An efficient numerical model has been developed and the related sub-physical models 
also have been incorporated for calculating transient two-phase combustion flows. 
The gas-phase equations are written in an Eulerian coordinate and solved using a 
control-volume based finite-difference method on an unsteady fashion. Spatial dif- 
ferences are formed on a curvilinear general coordinate with all gas phase variables 
stored at the same grid point. Second order accurate central differencing scheme is 
used for the diffusion terms and second order upwind scheme is used for the convec- 
tion terms. Whereas, the liquid-phase is presented in Lagrangian coordinates. This 
hybrid algorithm has flexibilities in handing poly-dispersed sprays in combustion 
devices and assigning individual particle properties. The two-way coupling between 
the two phases is described by the interaction source terms which represent the 
rates of momentum, mass and heat exchange. A strongly coupling numerical proce- 
dure based on PI SO algorithm with predictor/multi-corrector stages is developed 
and implemented into a computer code for multiphase all-speed transient flows in 
complex geometries (MAST). 
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A two-equation k — t turbulence model is used to characterize the gas phase 
mean and fluctuating flow properties. The particle turbulent modulation effect to 
gas phase is modeled within the stochastic framework. Several modulation models 
are tested and the numerical results are compared with available experimental data. 
Particle motion equation is the simplified B-B-0 equation and solved in Lagrangian 
coordinates for the present separated flow models. Particle dispersion by turbulent 
fluctuations is modelled by allowing particles to interact with a succession of eddies 
using a random-walk or Monte Carlo method. Parcel PDF model is used to im- 
prove computational efficiency. Poly-disperse is achieved by using size distribution 
function through Monte Carlo simulation or generated by breakup model. Droplet 
evaporation and heat transfer are calculated using Frossling correlation and Ranz- 
Marshall correlation respectively. Accurate calculation of mass and heat transfer is 
achieved by automatic reductions in the timestep when the exchange rate becomes 
large. The variable thermophysical properties are obtained from JANAF data bank. 
A turbulent diffusion flame and single step chemical reaction model is used for spray 
combustion. Dense spray effects are accounted for by a droplet breakup and collision 
model. Two breakup models including a Taylor analogy breakup (TAB) model and 
a wave instability model are incorporated and compared. The incorporated collision 
model stochastically calculates the outcome of every collision, either coalescence or 
grazing collision. 

Several particle dispersions in homogeneous and nearly-homogeneous turbulent 
flows are studied to calibrate the stochastic method for particle-turbulence inter- 
actions, and to evaluate the present parcel PDF model. The computations were 
performed for the solid particle dispersion in nearly-homogeneous turbulence and a 
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particle laden round jet in inhomogeneous turbulence. To account for dense spray 
effects, droplet breakup and collision models were employed for non-evaporating 
transient and steady sprays. Turbulence modulation effect was tested at a steady 
solid-cone spray. A transient hollow-cone spray showed complicated two-phase in- 
teractions. Finally, transient evaporating and burning sprays incorporating droplet 
evaporation and eddy-breakup combustion models demonstrated the sophisticated 
structures of such polydispersed sprays. 

All the calculations were performed on a CRAY X-MP/24 supercomputer. 

V.2 Conclusions 

The numerical models have been developed and tested for the analysis of dilute 
and dense spray- combustion flows. The major conclusions of the present numerical 
studies are as follows: 

1. A non-iterative numerical technique for computing time-dependent spray 
combustion flows has been developed and incorporated with existing so- 
phisticated sub-physical models. The method is a fully interacting combi- 
nation of Eulerian fluid and Lagrangian particle calculations. The inter- 
action calculations between the two phases are formulated on a pressure- 
velocity coupling procedure based on the operator-splitting technique. 
This procedure eliminates the global iterations required in conventional 
PSIC procedure and, hence, is efficient for transient calculations. 

2. A new set of constants for eddy life time and length scale have been 
evaluated. Comparing with dilute nearly-homogeneous particle turbulent 
dispersion experimental data, the present particle dispersion calculation 
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procedure has the flexibility of taking into account both of the gravity 
effect and the non-Stokesian drag law and gives more satisfactory result. 

3. For non-evaporating, evaporating, and burning dense spray cases, the pre- 
dictions show a reasonably good agreement with available experimental 
results in terms of spray penetration, drop sizes, and overall configu- 
ration of a burning-spray flame. However, quantitative differences exist 
especially at near nozzle exit locations. The discrepancies observed in the 
results are attributed mainly to uncertainties in the initial spray particle 
size and velocity distributions, the single-step fast chemistry employed by 
the combustion model, and the deficiencies of the k — t turbulence model 
dealing with the dense and combustion sprays. 

4. Three turbulence modulation models are evaluated and compared with 
standard k — e turbulence model and dense spray experimental data. The 
analysis indicates that the direct contribution of particles to gas phase 
turbulence properties is important for this dense spray case. This effect 
is attributed to the slip velocities between two phases at the fluctuation 
level. The addition of particles in the two-phase flows damps the turbu- 
lence intensities of the gets phase. Favorable agreement with experimental 
data is achieved when modulation effect is included properly. Predictions 
based on the model of Mostafa and Mongia [83] have a best agreement 
with the measurements of Wu et al. [120] . 

5. The present implementation of the parcel PDF model has successfully 
demonstrated the capability of accurately representing dispersion in nearly- 
homogeneous and inhomogeneous turbulent flows with improved efficiency 
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over the delta function SSF model. The numerical results also indicate 
that the present parcel PDF model has the capability of accurately repre- 
senting drop dispersion in dense spray situations with manageable number 
of computational parcels. 

V.3 Recommendations 

The present study has presented a comprehensive numerical model for spray com- 
bustion flows. There is, however, a wide scope for further improvements in both 
numerical aspects and physical models for spray combustion processes should be 
further investigated due to their complicated physical phenomena. The following 
recommendations are intended as suggestions for improvements and extensions of 
the present spray combustion modelling. 

1. Development of strong interphase coupling procedure by combining mul- 
tiple pressure correction procedure and Volume of Fluid (VOF) method. 
The present study does not include the effect due to the volume occu- 
pied by the particle phase which perhaps plays an important role in the 
atomization region. In liquid rocket engine case, the atomization process 
takes place over an extended region within the combustion chamber and 
sub-scale numerical models are needed to simulate this process. 

2. The droplet evaporation calculation procedure following KIVA-II [4] is 
used in this study. Droplets located in the same control volume are evap- 
orating one by one and see different intermediate gas temperature and 
vapor mass fraction. In high evaporating case, both gas temperature and 
vapor mass fraction change a lot in the calculation time step, that will 
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result in different evaporation history for each droplet even with same 
physical properties just due to different calculation sequences. This nu- 
merical procedure induced difference needs to be avoided through more 
precise numerical consideration. 

3. Incorporations of equilibrium and finite rate chemistry packages for effi- 
cient transient reacting flow calculations. The strong coupling between 
the flow fields and chemical reactions makes large variations in gas density, 
temperature and velocity etc.. More comprehensive operator-splitting 
techniques are required to accurately and efficiently predict reacting flows, 
especially for fast transient ones involving complex and stiff chemical ki- 
netics. 

4. More studies on turbulence modulation effect by particles, non-isotropic 
turbulence model such as the algebraic stress model and the second- 
moment closures for two-phase combustion flows are desirable. 

5. The present evaporation model does not include the supecritical temper- 
ature and pressure effect. In fact, this effect would be a major contributor 
to spray combustion dynamics in liquid rockets engines such as the space 
shuttle main engines (SSME) and advanced gas turbine combustors etc.. 
Droplets exposed to sufficiently high temperatures under supercritical 
pressures may be heated to the critical state or be gasified upon reaching 
critical point. Such mechanisms require extensive experimental inputs 
and the related researches are currently underway. At the same time, 
the present model does not take into account the group effect for droplet 
evaporation/combustion. In the group combustion model of Chiu et al. 
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[26], four combustion modes may occur for a droplet cloud. Bellan et al. 
[13] also proposed a group evaporation model to account for the dense 
spray effect. How to define the group parameters based on the spray 
dynamics and to incorporate these models into CFD code are receiving 
more researching attentions [54]. 



Appendix A 


Particle Turbulent Dispersion 
From a Circle 


P(x,y) = f exp(-[{z - Xo ) 2 + (y -yo) 2 ]/(2v 2 ) is)/ ids 

= iexp(-{(x- lo) 2 + (j/ - i/of]/(2<i J ) <fe}/(2irr 0 ) (A.l) 

Let’s change to cylindrical coordinate, 


X 0 = ToCOSOc 

yo = r 0 sina 
ds = roda 
x = rcosd 
y — rsinO 


(A.2) 


Substituting above relations into equation A.l yields 

1 f 2 * 

P{r,9) ~ 77 1? / exp{-[r 2 + - 2rr 0 cos(a - 0)\/(2cr 2 )}da 

[2,-Kcry J 0 


exp[-(r 2 + r 2 )/( 2a 2 )] 
(27rcr) 2 


r2TT 

I exp[rrocos(a — 9) / (2a 2 )]da (A. 3) 
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Substituting 0 = a — 9 into equation A. 3 yields 

P(rJ) = r~‘ e xp{ rr 0 cosP,SW 

(ZTay J_ e 

From residual theory of complex integration theory, we obtain 

/ 2rr—$ r ^ ^ 

exp(rr o cos0/<j 2 )dj3 = 2m f exp[rr 0 (z + 1/z)/(2ct 2 )] — 

■9 J \z \= 1 ^ 


= 27T 


/ [i 


rr Q 2 / rr 0 £ 

+ ^r + --- + ^f 
1! n! 


( rr Q z \ n rrg / rr Q \ n 

^4^- + • • -Hi + 44 + • • • + ^r- 


= 27r Ef 


V 2cr^ / 12 


(A.5) 


Hence, we obtain probability density function of particle position distribution for its 
turbulent dispersion from a circle with radius r 0 , 


p, „ _ exp[-(r^ + r;)/(2^)] f (3)' . 

' ’ 2™* 2J n! ’ 

n=0 


(A. 6) 


Note, this function is inependent of 9. For particle dispersion from a point source, 
where r 0 = 0, P(r, 9) is simplified to 


P(r,9) = 


exp[— r 2 /2cr 2 ] 
(27TCT 2 ) 


(A. 7) 
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